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ABSTRACT 

We present an extragalactic population model of the cosmic background light to interpret the rich 
high-quality survey data in the X-ray and IR bands. The model incorporates star-formation and 
supermassive black hole (SMBH) accretion in a co-evolution scenario to fit simultaneously 617 data 
points of number counts, redshift distributions and local luminosity functions (LFs) with 19 free 
parameters. The model has four main components, the total IR LF, the SMBH accretion energy 
fraction in the IR band, the star-formation SED and the unobscured SMBH SED extinguished with a 
HI column density distribution. As a result of the observational uncertainties about the star- formation 
and SMBH SEDs, we present several variants of the model. The best-fit reduced \ 2 reaches as small as 
2.7-2.9 of which a significant amount (>0.8) is contributed by cosmic variances or caveats associated 
with data. Compared to previous models, the unique result of this model is to constrain the SMBH 
energy fraction in the IR band that is found to increase with the IR luminosity but decrease with 
redshift up to z ~ 1.5; this result is separately verified using aromatic feature equivalent width data. 
The joint modelling of X-ray and mid-IR data allows for improved constraints on the obscured AGN, 
especially the Compton-thick AGN population. All variants of the model require that Compton-thick 
AGN fractions decrease with the SMBH luminosity but increase with redshift while the type-1 AGN 
fraction has the reverse trend. 

Subject headings: galaxies: active - galaxies: evolution - galaxies: luminosity function, mass function 



1. INTRODUCTION 

The extragalactic background light (EBL) represents 
the accumulated radiation generated over cosmic history 
after the Big Bang. The EBL in the X-ray, or cosmic 
X-ray background (CXB), is now known to be the relic 
emission of cosmic supermassive b lack hole (SMBH) ac- 
cretion (e.g. iComastri et aLlll995D . while the cosmic in- 
frared background (CIRB) arises mainly from the inte- 
grated emiss ion of cosmic obscur ed star formation ac- 
tivities (e.g. iCharv fc""lj lbaz 2001). Extragalactic deep 
surveys in the X-ray and IR bands have been key tools, 
identifying the discrete sources responsible for the CXB 
and CIB, respectively, and reconstructing the time-lines 
of SMBH accretion and dusty star formation. Cosmic 
energy generation in both X-ray and IR bands show a 
rapid increase back to redshift z ~ 1-2, with the density 
of lu minous sources generally peaking at a higher redshift 
fe.g.lCowie et alllToM ILe Floc'h et aT1l27)05TIUeda et al.l 
200.i IBarger et al]l2005j T This down-sizing character of 
evolution is thought to imply a mass-dependent forma- 
tion scenario for both SMBH and galaxies, with more 
massive systems formed at earlier times. This parallel 
behavior implies a statistical co-evolution of SMBH ac- 
cretion and star formation, with strong evidence for phys- 
ical co-evolution provided by the tight relationships be- 
tween th e masses of the stellar bulges and central BHs in 
galaxies (iKormendv fc Richstonelll995l : iMagorrian et al.l 
1199a iFerrarese fc Merrittil2000D . 

While rare in the local universe, IR luminous dusty 
sources are a major presence at high redshift (z > 0.7, 



e.g. ILe Floc'h et al.l [200l IMagnelli etHI l20ll . The 
quantification of fractional contributions by accretion 
and star formation to their luminosities is key to un- 
derstanding co-evolution and its relationship to mech- 
anisms, such as gas accretion or galaxy interactions 
and mergers, responsible for building todays galaxies 
and BHs. Studies of local luminous infrared galaxies 
(LIRGs; 10 11 L < L T ir < 10 12 L Q ) and ultra lumi- 
nous infrared galaxies (ULIRGs; 10 12 L Q < Ltir < 
10 13 Lq) have yielded major insights. Early on, op- 
tical spectroscopy established the importance of active 
galactic nuclei (AGN) in powering IR luminous galax- 
ies, and t he increase of this importance with the IR lu- 
minosity (IVeilleux et ail 119951 iSanders fc Mirabel 119961: 



lArmus et al.l 120071 IVeilleux et all 11999ft . X-ray surveys 
revealed these objects to be X-ray faint, and many 
of them are highly extinguished even in the hard X- 
ray band or completely C ompton-thick (jRisaliti et al.l 
120001: llwasawa et al.l 120 111) . Observations in the mid- 
IR with much lower extinction based on many avail- 
able AGN/star-formation diagnostics found that accre- 
tion contributes to the total luminosi ty about 10-20% in 
LIRGs and 20-50% i n ULIRGs (e.g . IGenzel et al.lll998l: 
Armus et all 120071: IVeilleux et all 120091 : iPetric et al.l 



20111 ). However, the situation beyond the local universe 
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is much more uncertain. Studies of sub-millimeter galax- 
ies (SMGs) at z ~ 2 show that the number fraction of 
objects hosting AGN or the SMBH luminosity fraction 
are lower than in local objects with similar l uminosities 
(I Alexander et al J 120051: iValiante et ail 120071; iPope et al.l 
2008). On the other hand, the mid-IR-selected hyper- 
luminous infrared g alaxies (HLIRG s) show significant 
AGN activity (e.g. IWu et all [20T1 . For LIRGs and 
ULIRGs, the selection method also plays an important 
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role in the derived SM BH luminosity fraction, ranging 
from about 5% to 30% (lYan et all 120071 IFu et al J 1201a 
iFadda et al.ll20loh . 

A complete census of the AGN population across a 
wide redshift range is necessary for a full understand- 
ing of co-evolution, but is hampered by the large col- 
umn densities often obscuring the AGN. Hard (> 2 
keV) X-rays arising from hot gas close to the cen- 
tral SMBH are efficient at identifying AGN, because 
the photon can penetrate a significant amount of gas 
(e.g. 70% photons at 5keV can get through 7Vhi=10 23 
cm ) to escape the galaxy. However, a poorly known 
number of Compton-thick AGN (defined as Nm > 
10 24 cm" 2 in this paper) escape detection because 
even hard X-rays are suppressed to undetectable lev- 
els. Current techniques to probe these objects com- 
bine X-ray data with ph otometry at other wavelengths 



20051 lAlonso-Herrero et al.l 120061 


Daddi et al.l 2007; 


Donlev et al.l l200St lAlcxandcr et al. 


20081 iFiore et al.l 


2009t iLuo et al.ll2011[). The dusty torus surrounding a 



SMBH absorbs the UV/optical radiation from the cen- 
tral accretion disk and the reprocessed radiation emerges 
in the IR and peaks at wavelengths significantly shorter 
than emission from star formation regions. The appar- 
ent gas-to-dust ratio of the AGN is significantly larger 
than the value in the Milky Way (e.g . lRisaliti et alj|200a 
iMaiolino et alJl200HlShi et al.ll2006l ). leading to low mid- 
IR optical depths even for Compton-thick AGN. Still, 
broad-band mid-IR AGN selection techniques suffer from 
both contamination by star-forming galaxies and incom- 
pleteness in AGN identification (e.g. IDonlev et al.ll2008t 
iPetric et alJl20H . 

Galaxy population models have been used to represent 
information in deep surveys and EBL, typically focus- 
ing on either AGN a nd CXB Ce.g. ISetti fc WoltierlllQSl 



Madau et all 11994 IComastri et al.l 119951: lUeda et all 



20031: iTreister fc Urrvl 120051: IGilli et"alTl200lL I2007T) 



1991: Chary & Elbaz 2001: Xu et al. 2001: Lagachc et al. 


20031 12004: iRowan- 
20091: fVaiiante et al. 


Robinson 2009|: iLe Borgne et al.l 
2009: iFranceschini et all 12010b 


Gruppioni et al.ll2011 


;lBethermin et al.ll2011l). However, 



the rich collection of multi-wavelength data accumulated 
during the last decade make coherent models taking into 
account simultaneously star formation and SMBH ac- 
cretion not only possible, but necessary. This need is 
demonstrated by several lines of evidence that Spitzer 24 
fim sources have significant AGN contributions to their 
mid-infrared luminosity. At the high flux end (/24/nm > 5 
mJy), the Spitzer 5MUSES legacy program (PI: George 
Helou) , a mid-infrared spectroscopic survey of 330 galax- 
ies selected essentially by their 24 /xm flux, found 30% of 
its spectra have a low 6.2 /xm ar omatic equivalen t widt h 
indicative of AG N dominance dWu et al.l IIM |20TTI ). 
At lower fluxes, IDonlev et al.l (|2008l ) employed multi- 
ple broad-band photometric methods to look for sources 
with AGN-dominated mid-IR emission, and found that 
the fraction varies from 30% at fi^m — 1 mJy to around 
10% at /24/im=0.2 mJy and remains at 10% down to the 
survey limit (/24 M m=0.08 mJy). These fractions are gen- 
erally consistent with the Sp itzer spectroscopic results 
of relatively small samples (|Yan et al.l 120071 : IFu et al] 



I20iq IFadda et al]|2010D . Further demonstration of this 
need comes from Chandra surveys, currently reaching 
deep enough to explore X-ray emission from star forma- 
tion. The fraction of galaxies powered by star formation 
rises steeply with decreasing X-ray flux, reaching around 
50% at the survey limit (10" 17 erg/s/cm 2 at 0.5-2 keV) 
(|Bauer et al.l 12004). Since these fainter fluxes are also 
populated by heavily obscured AGN, taking into account 
the star formation contribution is critical for deriving a 
correct AGN census. Attempts to combine X-ray and IR 
data to improve the constraints on th e obscured AGN 
have produced some in teresting results (Ballantvn e et al.l 
120061 : lHan et al.l [20121 ). indicating higher obscured frac- 
tions at higher redshift. 

This paper introduces a model of galaxy populations 
across cosmic time that fits simultaneously X-ray and 
IR survey data and reproduces the X-ray and IR extra- 
galactic background. This approach allows estimations 
of quantities inaccessible with either IR or X-ray data 
alone, such as the SMBH energy fraction as a function 
of the IR luminosity out to z ~ 3. In this first paper, 
we present the models construction and the basic out- 
puts, namely the total IR Luminosity Function (LF), 
the distributions of SMBH energy fraction and HI col- 
umn density. Our philosophy is to incorporate known 
observational results in the model instead of leaving ev- 
erything as free parameters, but also present several dif- 
ferent variants to reflect the uncertainties on the adopted 
observational results. In addition to this paper that fo- 
cuses on the model construction, we will discuss in other 
two papers the implications on the Compton-thick AGN 
abundance and the cosmic link between obscured star 
formation and SMBH accretion, respectively. The overall 
structure of this paper is summarized below. We briefly 
summarize previous models in § [51 In § [31 we describe 
the structure of our model. The general procedure to 
fit the observed data is given in § [U The datasets used 
for the fit are listed in § [51 The results of the fit and 
the best-fit parameters are shown in § [6j The fundamen- 
tal outputs of the model are discussed in § [TJ Through 
out the paper, we divide galaxies into three types based 
on the X-ray luminosities: quasi-stellar objects (QSOs) 
with logL 2 -iokcV > 10 44 erg/s, Seyfert galaxies with 10 42 
erg/s < logL 2 -iokcV < 10 44 erg/s, and BH quiescent 
galaxies (BHQGs) with logL 2 _iokcV < 10 42 erg/s. AGN 
include QSOs and Seyferts. Based on the HI column den- 
sity, the AGN are divided into type-1 (7V H i < 10 22 cm -2 ); 
Compton-thin type-2 (10 22 cm" 2 < N m < 10 24 cm" 2 ) 
and Compton-thick {{N m > 10 24 cm -2 ). The total IR 
luminosity is defined to be from 8 to 1000 /im. The galax- 
ies are also defined according to the total IR luminosity as 
IR quiescent galaxies (IRQGs; Ltir < 10 11 L Q ), LIRGs 
(10 11 L < Ltir < 10 12 L ), ULIRGs (10 12 L Q < L Tm 
< 10 13 L ) and HLIRGs (Ltir > 10 13 L ). We adopt 
a cosmology with Hq=70 km s" 1 Mpc" 1 , £! m =0.3 and 
Q A =0.7. 

2. BACKGROUND: THE SYNTHESIS MODEL OF 
EXTRA-GALACTIC BACKGROUND LIGHT 

Before we describe our EBL model, we here summarize 
previous CXB and CIRB models. In general, a synthesis 
model first parametrizes the LF at a wavelength and its 
redshift evolution. With assumed SEDs, the model can 
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then predict a LF at any observed monochromatic wave- 
length, which will be compared to the data to derive the 
best-fit parameters. The approach can, in principle, fit 
all the survey data regardless of field size, survey depth 
and wavelength coverage. 

2.1. Cosmic X-ray Synthesis Model 

The CXB was dis covered almost half a century ago 
(jGiacconi et al.lfl962t ) , and is now understood to be dom- 
inated by the integrated emission from gas being heated 
and accreted onto SMBHs at the centers of active galax- 
ies. The pure black-body spectrum of the CMB places 
strong upper limits on the contributions to the CXB 
from diffuse hot intergalactic gas (10 -4 , iWright et al.l 
1994). Deep extra-galactic surveys have now resolved 
a significant fraction (~70%) of the ba ckground light 
at < 10 keV into discrete s ources (IHasinger et all 
19981: iMushotzkv et all 120001: lAlexander et all 120031: 
Worslev et al.l I2005D . with the resolved fraction de- 



creasing with the increasing energy bands. The exact 
fraction also depends on the absolute normalization of 
the CXB spectrum whose dispersion reaches ~30% at 
low energies (< 10 ke V) and ~10% above 10 keV among 
differ en t studies (e.g . IMarshall et alj|l980j: IGruber et all 



19991: iVecchi et all 119991: iReynivtsev et al.l 120031: 
Aiello et al.l 120081: iChurazov et al.l l2007t iMoretti et al l 
■-'1)09: iTurler eTalll2010l) . ~ 



The CXB spectrum peaks around 20-30 keV and drops 
toward both lower and higher energies, and its slope at < 
20 keV is harder than those of Seyfert-1 galaxies. In order 
to reproduce this slope and turnover around 20-30 keV, 
the general strategy is to invoke the SED of a type-1 AGN 
along with a varying HI column density (iSetti fc Woltierl 
U989HMadau et all 119941: IComastri et all 119951). consis- 
tent with the unified scheme for AGN ( Antonucc11ll993l) 
where Seyfert-2 galaxies are obscured variants of Seyfert- 
1 galaxies and exhibit harder X-ray spectra. Today 
CXB models not only reproduce the CXB spectrum (1- 
100 keV) but also the X-ray number counts, redshift 
distributions, an d spectral properties in both soft and 
hard bands (e.g. ITreister fc Urry| |200"5l: IGilli et al.l 120071: 
iTreister et all I2009D . It is generally accepted that the 
majority of X-ray AGN activity occurs in the obscured 
phase and that the evolution of X-ray AGN is luminosity 
dependent, with the density of Seyfert-like AGN peak- 
ing around z=l and quasars peaking at z=2-3. How- 
ever, there are still important uncertainties in the de- 
rived parameters, such as the type-l/type-2 AGN ra- 
tio and its variation with luminosity and redshift, and 
the amount of Compton-thick objects. Such difficul- 
ties are at least partly caused by the fact that the cur- 
rent X-ray surveys at < 10 keV are not sensitive to 
heavily-extinguished AGN, especially Compton-thick ob- 
jects, and overall number counts remain limited. 

2.2. Cosmic Infrared Synthesis Model 

The CIRB was first dire ctly measured in the 19 90s by 
COBE (for a review, see lHauser fc Dwekl l200l . The 
energy of CIRB is comp arable to the cosmic UV/optical 
energy (|Dole et al.ll2006 ) , while the ratio in the local uni- 
verse is only one third ( Driver et al.ll2008h . This implies 
an increase in dusty sources with increasing redshift. The 
slope of the CIRB in the far-IR/submm range is flatter 



than that of local starburst galaxies, which provides fur- 
ther evidence for the strong evolution of IR sources whose 
peak emission is redshifted toward longer wavelength, 
flattening the CIRB spectrum at the longer wavelength. 

Unlike the CXB model where the HI column density 
is considered to be the main factor driving the varia- 
tion of the SED, the CIRB model typically varies the 
SED mainly as a function of IR luminosi ty. They either 
adopt a fixed SED at a give n luminosity ()Charv fc Elbazl 
l2001tlLe Borgne et all20 09). or further c onsider the scat- 
ter of the SED at a give n luminosity (jChapman et all 
120031 : iValiante et all 120091 ) . or assume different IR pop- 
ulations that exhibit dif ferent SEDs and occupy dif- 
feren t luminosity ranges dXu et all 120011: lLagache et al 
20031 12004 iRowan-Robinsonl 120091: iFranceschini et al 



2010: Gruppioni et al. 2011; Bcthcrmin et al. 2011). The 



SMBH radiation, with typically a much warmer IR SED 
compared to that of star formation regions is energeti- 
cally important in mid-IR surveys: the contribution of 
AGN to the mid-IR background reaches noti ceable frac- 
tions, e.g. 10-20% in the 1-20 u m. range by iSilya et all 
p00l or 2-10% in the 3-24 /mi bv ITreister et all (|2006f) : 
classifications (AGN vs star-formation dominated) of 
24 /im sources indicate that the fraction of objects 
whose mid-IR emission is dominated by AGN ranges 
from ~40% at f24 Aim > 5 mJy and drops down with 
decreasing fluxes but remains aro und 10% down to 
the faintest end f24/;,m ~ 0.08 mJv (iDonley et ail 120081 : 
IFu et al.ll2010l : IWu et all I20T0I: IChoi et all 12011^ How- 
ever, the 24 /im number counts can be reproduced 
equally well by models w it h AGN (IValiante et all 120091: 



IFranceschini et aTl 120101: iGruppioni et al.l 1201 ID and 



witho ut AGN ([Bethermin et all 



Riom et a 
20111: ID 



c Borgne et al. 



2009), reflecting the limitation of the CIRB synthesis 
models in decomposing SMBH and star-formation emis- 
sion. Models without AGN can re-produce the effect 
of AGN in the mid-IR simply by modifying the star- 
formation mid-IR SED or by adjusting the slope of LFs 
and their redshift evolution. 

3. THE CONSTRUCTION OF A JOINT MODEL OF X-RAY 
AND IR BACKGROUNDS 

Our model has four basic components, each of which 
contains a few components. The four components are a 
luminosity function defined in the total IR band (8-1000 
/mi), the SMBH energy fraction in the total IR band, the 
luminosity-dependent SED for the star-forming compo- 
nent, and the type-1 SMBH SED extinguished by a HI 
column density. The first component describes the num- 
ber density of galaxies at a given redshift and total IR 
luminosity, while the second one decomposes this lumi- 
nosity into star-formation and SMBH total IR emission. 
The last two then generate star-formation and SMBH 
emission at any observed band from their corresponding 
total IR luminosities. The general rule for setting the 
total number of free parameters is that the model starts 
with the minimum set of free parameters and additional 
ones are added if the reduced x 2 decreases by about 0.5. 
The goal is to fit model predictions to number counts, 
redshift distributions and local LFs data, which offers 
hundreds of data points. As discussed in the following 
of this section, we will run a total of four variants of the 
model to reflect current uncertainties on the star-forming 
and SMBH SEDs. 
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Fig. 1. — The distribution of the SMBH energy fraction in the 
total IR band of 74 1 -Jy ULIRGs based on the aromatic feature 
EW measurements by Veillcux et al. (2009). The filled histogram 
shows the derived SMBH energy fraction distribution. The solid 
line is the Gaussian fit to the filled histogram. For the detail, see 



3.1. Infrared Luminosity Function 

We represent the galaxy total IR LF as a double power- 
law distribution with ten free parameters: 

\ dN 



dlogL x dV 



(£/£*) 71 + ( L / L * 



[Mpc^logLQ 1 ] (1) 



We adopt a simple independent luminosity and density 
evolution (ILDE). The evolution is described as a poly- 
nominal function of redshift: 

logi*^) = logL^o + V x h,l + V 2xk 2,i + V 3x h,i (2) 

log$*(z) = log^o +T)xki td + ?7 2 xfc M + ?7 3 xfc 3 ^ (3) 

where r\ — log(l + z). For this component of the model, 
there are ten free parameters including four in the lumi- 
nosity evolution (logL^o, &i,z, ^2,z & ^3,i)j f° ur m the 
density evolution (log^o, &2,d & ^3,d) and two 

slopes (71 Sz 72). 

3.2. The SMBH Energy Fraction 

Individual objects in the model are assumed to be ex- 
periencing both star formation and SMBH accretion ac- 
tivities, so that Ltir = ^tir + ^tir- F° r the formula of 
the SMBH energy fraction, we assume a bounded Gaus- 
sian distribution for the logarithm of the SMBH energy 

fraction in the total IR band / B : B = ^™ , as motivated 
by works of iValiante et al.1 (|2009f ): 



BH 



dN 



cflog/, 



BH «exp(- 

TIR 



(log/ BB -log/ BB R ) 2 . 



2(ct bh ) 2 



P^log/^ R )dlog/^ R 



(4) 
(5) 



where the boundary of log/rj?^ is set to have the min- 
imum and maximum at -6.0 and 0.0, respectively (for 
details, see § [4j) . 

This Gaussian distribution is assumed to be redshift 
and luminosity dependent. The mean value of the distri- 
bution (log/? ) is a function of both total IR luminosity 
and redshift: 



log/? H = log/ 



BH 



- Pf E 



where 



log/™ = iog/,7 + v xkZ3 + fxkft 



(■BH 



: xloe 



BH 



^TIR 



■ BH 



(6) 



(7) 



The width of the distribution is only luminosity depen- 
dent: 

cr BH = ( T BH +^xlogi^L (8) 

At z—0 and Ltir=£*,/bh, the SMBH energy fraction 
is fixed at / BH , which is based on the sample of local 
1 Jy ULIRGs' (|Kim fc Sandersl[T99l . The advantage of 
ULIRGs compared to lower luminosity objects is that 
the SMBH contributions are important so that the dis- 
tribution of the SMBH energy fraction can be derived 
to a high accuracy. For a sample of 74 ULIRGs, we de- 
rived the SMBH IR luminosity from the observed 6.2 
/zm equivalent width and 6.0 / zm co ntinuum luminosity 
as measured bv iVeilleux et al.l fj2009f ) : 

L e!L = max[(1.0 - EWe.2^m/EW Bt3 , T -f OTina ,ti n) x 

L6p m xexp(T 6lim ),0] (9) 

where EWe.2fj, m is the observed 6.2 /im equivalent width, 
Ls^m is the observed 6.0 /zm continuum luminosity, 
-EWstar-formation is the equivalent width of pure star for - 
mation and taken to be 4.2 /an (|Veilleux et alJl2009T) . 
and Te^m is the effective optical depth at 6 /im and set 
to be 10% of the measured 9.7 um s ilicate depth using 
the SMC extinction curve (jPel Il992h . Given the star- 
forming EWs still have a small scatter (e.g. iWu et ahl 
120101: iStierwalt et~alll2012t) . we set the maximum value 
of the derived SMBH radiation in Equation[3]to be zero. 
Note that the aromatic feature EW relies on the method 
to measure it, e.g., sometimes the wing emission of the 
feature is included through extrapolating the analytic 
profile. As long as the same method is adopted, the 
SMBH energy fraction from the above equation should 
not differ significantly. The is then corrected to 

the total SMBH IR luminosity by a factor of 1.6 using 
our own SMBH SED. Combined with the measured total 
IR luminosity, the final distribution of the SMBH energy 
fraction is shown in Figure [T] Our assumption of a Gaus- 
sian distribution is more or less correct given the derived 
shape. The fit with a Gaussian gives log/?^=-1.21. The 
median luminosity of this sample gives logL* jbh=12.3. 
Note that the above values do not depend significantly on 
the adopted SW star -formation since EWe^m is small in 
the majority of the objects. However, we did not use the 
above distribution to fix tr BH but set it as a free parame- 
ter, because the width measured in Figure[T]is broadened 
by the intrinsic scatter in EWsf and luminosity correc- 
tion Ltir/ Ls^m. We notice that the median fraction de- 
rived here is significantly lower than the energy fraction 
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of roue hly 35-40% as derived by iVeilleux et all ((2009() . 
The main reason is that their numbers are given in the 
bolometric band while ours are integrated from 8-1000 
fim. While the emission is similar between these two 
bands for ULIRGs, a factor of ~ 5 needs to be applied 
for the SMBH radiation. 

There are five free parameters in total, including two 
describing redshift evolution of the mean fraction (kf^ 

& fcf d)) one characterizing luminosity dependence of the 
mean fraction (p^bh), one for the scatter (cr^ H ) at the 
total IR luminosity of 10 12 3 L Q and another for its lu- 
minosity dependence {p a )- 

3.3. Star Formation SED Template 

For the IR SED of the star-formation emission, we 
adopted the luminosity-dependent SED shapes and as- 
sociated dispersions. This set of SEDs is well character- 
ized in the local universe; the uncertainty lies at high 
z where the difference is already observed but not fully 
quantified. 

3.3.1. Local Infrared SED Templates 

It is now established that the IR luminosity is not the 
main driver of the IR SED shape, i.e., at a given luminos- 
ity, t he IR color shows a lar ge scatter (|Chapman et al.1 
2001 IHwang et al.1 120101 : iRujopakarn et al.l l20liaF 
Elbaz et al.1 201 lK Fortunately, this scatter has been 



well quantified. We thus adopted the observed trend be- 
tween IR luminosity and color ( /^mAif- foo^'m) ) as we ^ as 
the as sociated scatter as established by [Chapman et ahl 
( 2003) for local IRAS galaxies. The probability distribu- 
tion of IR color at a given star-forming IR luminosity is 
described by a Gaussian function: 



P 



SF-color 



dN 

dC 



exp(- 



(C - Cq) 2 
2(a c ) 2 



where C=log( j- 



^«), a c =0.065, and 



C =C, +7xlog(l + 



L 



TIR * 



5xlog(l 



L 



*,c 



L 



(10) 



) (11) 



TIR 



C* = -0.35, 7=0.16, 5=0.02, L^c=5.0xl0 10 L , and 
L TIR is the star-forming total IR luminosity from 3 /ini 

to 1100 ^m. For a given color of /"//^s - loo am) ' we 
used t he infrared SED template library of iDale & Heloul 
(|2002t hereafter DH02), which sorts the IR SED by 
that color. The template was built empirically first 
from IRAS/ISO 3-100 /im observations (jDale et al.ll2001| ) 
and then extended up to submm wavelen gths based on 
SCUBA and ISO data (jDale fc Heloull200l . Spitzer data 
(3-160/im) confirms its ability to cover the overall range 
of SED variations and describe the med ian trend between 
different colors below 200 [im (e.g. IDale et all 120071: 
IWu et ahl I2010D . Recent Hersc hel far-IR observations 
furth e r demonstrated its p ower (I Rowan- Robinson et ahl 
l2010t IBoselli et al.l I20TI IDale et all 120121) . Never- 
theless, observations suggest that this set of tem- 
plates likely lacks the dynamic range in the SED 
variations above rest-frame 200 /xm. For example, 
DH02 templates do not produce objects with rest-frame 



f v (Herschel-500jnm)/,f 1/ (Spitzer- 24^m) > 6.5 as observed 
in iRowan-Robinson et al.l (|2010f ) although they bracket 
well the spread in /^(Spitzer-70/im)//,,(Spitzer-24^m). 
The templates also under-repre sent the diversity of slopes 
in the mid-IR (|Wu et al. 2010j). It should be emphasized 
that all template libraries in the literature are based on 
broad-band photometry, mainly using ISO and IRAS fil- 
ters. However, the detailed SED shape at fine scales 
within a filter can actually affect model predictions for 
broad-band surveys as the different parts of spectra are 
rcdshifted into a fixed wavelength filter. 

3.3.2. IR SEDs At High-z 

The sensitive infrared surveys carried out by 
Spitzer and Herschel are revealing the SED shapes 
of high redshift luminous sources to be signifi- 
cantly differen t from local galaxie s at the same lu- 
minosity (e.g. iPapovich et all I2IM iRigby et al.l 120081: 
I Svmeonidis et al.ll2009l: IMuzzin et al ll2010UHwang et al.l 
l2010t IRujopakarn et alj l2011aUbl : lElbaz et al.l 120111) . 
Many high-z LIRGs and ULIRGs show IR SEDs similar 
to those of local less luminous galaxies, i.e. high-z LIRGs 
and ULIRGs are colder compared to their local counter- 
parts; IR luminosities do not correspond to unique prop- 
erties of star-forming regions. To account for this fact, 
we evolve the L* c m Equation 1111 through 



L* tC (z)=L*, c (z = 0)(l + z) p ° 



(12) 



It is, however, still difficult to constrain p c because of 
the limited number of available broad-band photometric 
bands, the dynamic range of the IR luminosities at a fixed 
redshift and the observational selection biases. It is pos- 
sible that the evolution in the star-formati o n SED only 
takes place at z > 1.5 (Hwang et a.1.1 f20Tol lElbaz et al.l 
I2010T ) . But it is still not certain whether the basic shape 
of Equation I 111 i.e., a two-power-law mean trend and 
a Gaussian distr ibution of the scat ter, is valid at high- 
z. For example, lElbaz et al.1 (|2011[) have described the 
high-z LTm/Ls^m color as a Gaussian distribution plus a 
higher- value tail that is independent of the IR luminosity. 
In this paper, we assumed that the form of the local rela- 
tionship holds at the high z and explored two possibilities 
of p c , zero for no evolution and 2.0 for rapid evolution 
to bracket the range as discussed in the l i terature (e.g. 
I Svmeonidis et alj|2009t IMuzzin et al.ll2010l: IHwang et al.l 
hoiOt iMagnelli et al.l l2009t iChapin et al.ll2011l ). Notice 
that we did ran the model with various p c values from 
0.0 to 3.0 but failed to found that a non-zero p c value 
can give the better fit compared to the no SED evolu- 
tion case. Because of this, the p c =2.0 is used to illustrate 
the effect of the SED evolution in the model's fitting. 

3.3.3. X-ray SED 

Deep X-ray surveys, especially in the soft band, can 
probe emission from star- formation regions. We thus re- 
lated the star-forming X-ray emission to the star forma- 
tion IR luminosity and constructed the X-ray spectrum 
to go with each IR SED template. In galaxies with active 
star formation, the X-ray emission is dominated by short- 
lived high mass X-ray bi nary (HMXB) and is known to be 
relate d to the SFR (e.g. iGrimm et al.l l2003l: iPersic et al.1 
120041 : iLehmer et al.ll2010f h In quiescent objects, the low- 
mass X-ray binary (LMXB) emission becomes important 



6 



Shi et al. 





Fig. 2. — The relationship between IR and X-ray emis s ion o f 
star-forming galax ies ba sed on the samj 3le of IPersic efaTI 1(2001 ) 
and ILeh mcr ct al. (2010) (see Equation 1131 ■ 



and X-ray luminosity is no longer a reliable SFR tracer. 
Nevertheless, IR emission also has a contribution from 
dust heated by long-lived stars, so it is reasonable to 
expect a relation between IR and X-ray even in qui- 
escent objects. W e th us compiled measur ements from 
IPersic et all (|2004t l and lLehmer et al.l (|2010T > but exclud- 
ing AGN and AGN+SF composite objects. For a few 
common objects in both studies, data from the latter are 
used. The total IR luminosities unavailable in the former 
study are measured by our own based on the IRA S four 
band fluxes through the lSanders fc Mirabell (|1996D equa- 
tion. The final sample covers the IR luminosity from 10 s 
up to > 10 12 L Q , as shown in Figure [2] A two power law 
fit to the data gives: 



wovelength(//m) 



Fig. 3. — The IR spectrum of type-1 AGN based on the Spitzer 
IRS spectrum and MIPS photometry of PG quasars (the upper 
bound of the grey area). First, the median IR SED is measured 
for the whole PG sample. The star formation contribution is sub- 
tracted based on the measured aromatic features of this median 
spectrum. The black solid line and filled circles are the IRS spec- 
trum and MIPS 70 fim photometry after subtractions. Note that 
the observed MIPS 160 fim is dominated by star formation. The 
purple dotted line describes the underlying continuum of the sil- 
icate features. It is a fit to the 5-8fim IRS spectrum and 70 fim 
photometry with two power laws where the index at long wave- 
length is fixed to be 4 to mimic the modified black-body emissivity 
of 1.0 [yfu oc u 4 ). A single power law (red dashed line) connecting 
the IRS 30 fim emission and MIPS 70 /an photometry is used to 
model the silicate emission above the continuum roughly between 
30 fim and 38 (mi. 

The grand unified model of AGN (jAntonuccil 119931) 
states that all AGN are intrinsically the same and the 
emerging SED is only determined by the line of sight with 
respect to the dusty torus; the type-2 objects are the ob- 
scured version of the type-1 objects. Adopting this sce- 
nario, we first constructed the intrinsic AGN SED from 
type-1 data-set and then modified it with extinction as 
characterized by the HI column density (Njji) to obtain 
the extinction-dependent SED. A reflected component in 
the X-ray is always added: 



lo g-^2-10kcV 



-4.70+ 1.08 x logi|f R 
-0.59 + 0.68 x logL|f R 



for logipa < io.3l 



for log- 



> 



exp(- 



f 

-t(v, N m )) 



AGN-obscrvcd 



pAGN-intrinsic . 



Ju 

reflected — component 



(14) 



l q 
(13) 

where all luminosities are in units of solar luminosity. 
The 1-tr scatters for trends at the low and high lumi- 
nosity ends are 0.28 and 0.27 dex, respectively. The 
fit does not use the only one upperlimit but the re- 
sult does not chang e if including it. As discussed by 
iLehmer et al.l (|2010D . the scatter between X-ray emis- 
sion and SFR can be reduced if including the stellar 
mass quantity which, however, is not feasible for our 
model that does not have optical/near- IR components 
to infer stellar masses. Independent of the X-ray lumi- 
nosity, the spectrum from 0.2 to 100 keV is fixed to be 
that of the star-forming galaxy M82 based o n the Bep- 
poSAX o bservation from Cappi et al.l (|1999| ) and Swift 
data fromfBaumga rtner et al. ( 20111 ). We further assume 
no redshift evolution for the ^x^ray — relationship 
and X-ray spectrum. The lack of redshift evolution in 
the -^x-ray — ^iiC relationship is tentatively suggested by 
some works (| Grimm "eTaT1l2lMrSvmeonidis et al.ll20ll . 

3.4. SMBH SED Template 



10.31 

where /AGN-obscrvcd ig the orjserv ed SED of AGN, 
jAGN-mtnnsic - g ^ ne intrinsic SED before extinction and 
t(v, Nm) is the extinction curve. 

The significant difference from the star-formation SED 
templates is that we only used the mean SMBH SED but 
did not consider its dispersion. This is partly because to 
keep computations manageable but also because of the 
lack of studies that address fully the dispersion of the 
AGN SED from X-ray all t he way to submm as a function 
of the BH luminosity (e.g. IHao et al.|[2012n . 

3.4.1. The IR/submm SED of Unobscured AGN 

The distinct IR continuum feature of the unobscured 
AGN is strong emission of dust heated to near the subli- 
mation temperature (~ 1000 K) by UV/optical radiation 
from the accretion disk. This feature serves as the basis 
of va r ious mid-IR AGN selection techniqu e s (lLacv et al 
2004 1 : IStern et alj|2005t lArmus et al.l 120071 : iDonlev et al 
20081 IWu et all 1201 ID . The shape of this feature has 
been well characterized with the emission peaking at 3-5 
/im and the 1-30 /im slope a=l±0.3 (/„ oc v~ a \ al- 
though the slope estimate depends on the definition of 
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Fig. 4. — The relationship between 2-10 keV unobscured X-ray 
and 12.3 /im luminosities. The dashed line is the best fit for all 
Seyfert and LINER galaxi es with well- resolved images at 12.3 fim 
(circled symbols; Gandhi ct al. 2009). The solid line is the fit 
adopted in this paper that is exclusively for PG quasars. It is 
derived by fixing the slope as that for Seyfert and LINER ones but 
recalculating the normalization to produce the median trend of PG 
quasars. The offset between the dashed and solid lines is due to 
the strong silicate features in the PG quasars. 

the photometric bands and their placement with respect 
to the silicate feature besides the intrinsic dispersion. 
Generally, the broad-band photometr ic observations in- 
dicat e slopes between 1 .0 and 1.5 (e.g. iNeugebauer et al.l 
119871 : lHaas et all 12003). The Spitzer IRS observations 
of a large sample of local AGN reveal smaller slopes 
for the continuum underlying the silicate feature. For 
example, the whole sample of 87 PG quasars shows a 
relatively universal slope with median and sta ndard de- 
viations of a5 -30um=0-75± 0.35 ( f v ex u~ a ) (IShi et al.1 
I2007T ) (also see lNetzer eFaT1l200rt [Maiolino et al.ll2007l ). 
The Seyfert 1 objects i n the complete 12 (im sample have 
ai5-30 M m=0.85±0.61 (|Wu et al.ll2009l ) where the slightly 
steeper slope and the larger dispersion are likely caused 
by the larger star formation contamination in these low 
luminosity AGN. No obvious evidence has been seen for 
the dependence of the slope on the BH luminosity in the 
above studies. 

The most prominent spectral feature of the unobscured 
AGN SED is the silicate emission features at 9.7 and 18 
/im as established unambiguously by Spitzer observations 
(e.g. lHao etaHl200alSni et al.l l2006) . It is important to 
include this feature in the SED template for our model as 
it can affect the broad-band photometry The strength 
of this feature depend s weakly on the SMBH luminosity 
(|Maiolino et al.ll2007|) . with only a factor of 2 increase 
for a SMBH luminosity range of 10 4 . We thus adopt the 
same feature strength independent of the SMBH lumi- 
nosity. 

The most uncertain part of the SMBH IR SED lies 
above ~ 30 /im where the star-formation emission be- 
comes important or even dominates for quasars. One so- 
lution is to spatially separate the nuclear emission from 
star formation regions in host galaxi es; such obseryation s 
are only feasible below 30 /im (e.g. IGandhi et aLll2009h . 
Another challenging problem with this strategy is that 



we have no idea if the dusty torus itself is experiencing 
star formation and how important that is. The current 
approach is to scale the star formation with the observed 
aromatic feature and the n subtract it off (e.g. IShi et al"] 
I2007t INetzer et all 120071) . This leads to the conclusion 
that above ~ 100 /im, all the observed emission typically 
comes from star formation, and that SED part then re- 
lies on the extrapolation from the short wavelength. In 
spite of significant caveats associated with the AGN SED 
above 30 /im, its effect on our model is limited. This is 
because the majority (90%) of the AGN IR emission (3- 
1000/xm) emerges below 30 /xm; their far-IR contribution 
is negligible to the observed far-IR radiation. 

In practice, we constructed the unobscured IR SED 
based on Spitzer observations of the whole PG sample 
that is selected based on the UV e xcess and are thu s 
type-1 objects with little obscuration (IGreen et al.I 1986). 
To construct the intrinsic SED, we first derived the me- 
dian composite spectrum of Spitzer IRS spectra (5-30 
(im) and MIPS photometry (24, 70 and 160 jum) of all 
the PG objects. With the composite spectrum, we mea- 
sured the 11.3 /zm aromatic feature and scaled it to de- 
rive the star-forming emission at the Spitzer wavelength 
based on the DH02 template. The star-formation emis- 
sion contributes less than 15% at < 30 /im and reaches 
^50% in the 70 /mi filter while dominating the 160 /im 
filter emission. A continuous IR SED template is then 
constructed based on the star-formation subtracted IRS 
spectrum and photometry in the 70 /im filter where all 
the aromatic features are also carefully removed. To fill 
the gap in the SED beyond 30 fan, as shown in Figure [31 
the underlying continuum of the silicate feature is fitted 
with a double power law (purple dotted line) to the 5- 
8 /tm IRS part and MIPS 70/tm photometry where the 
power law index at long wavelength part is fixed to be 
4 to mimic the modified black-body with emissivity of 
1.0 (vf u oc i/ 4 ). A single power law connecting the IRS 
30/xm emission and MIPS 70/im photometry (red dashed 
line) is used to model the silicate emission between 30 
and 38 /tm. The resulting IR SED is quite consistent 
with the p hotometry based SED derived from the SDSS 
quasars by IRichards et all (|2006| ). except for the much 
better characterized silicate features in our template. 

Extracting the intrinsic IR SED for a large sample of 
AGN with a range of host galaxy contamination, extinc- 
tions, etc, would be quite difficult. Moreover, the spatial 
and spectral resolution data are not available to allow 
for a result comparable to that for star-forming galaxies 
with well characterized mean SED and associated scat- 
ter as a function of th e luminosity ([Dale fc Heloul [20021 
IChapman et al.l l2003f) . We thus adopt the mean PG 
quasar IR SED with no luminosity dependen ce or scat- 
ter. The recent work by Mullanc v et al.l (|2011D claimed a 
luminosity dependence of the AGN intrinsic mid-IR/far- 
IR ratio in a small AGN sample. However, their low- 
and high-luminosity sub-samples have completely differ- 
ent extinctions, which could easily induce such a depen- 
dence. The studies of X-ray AGN in the COSMOS field 
show no ap parent luminosit y and redshift dependence for 
their SED i|Hao et al.|[20ll . 

3.4.2. X-ray SED As A Function Of N m 

For the X-ray SED as a function of the HI column den- 
sity, we basically adopted the work of lGilli et al.l (2007). 
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Since no scatter is included for the IR part in our model, 
we did not consider the scatte r in th e X-ray spectral 
slope either. While iGilli et al.l (|2007l ) argued that the 
scatter in the X-ray spectra is important to constrain 
the Compton-thick AGN population, this will not be the 
case for our model due to the fundamental difference in 
the way that the Compton-thick AGN abundance is con- 
strained. Their CXB model derived the Compton-thick 
AGN by subtracting the Compton-thin AGN contribu- 
tion to the CXB spectrum around 30 keV. The extrap- 
olation of Compton-thin AGN contribution from < 10 
keV to 30 keV certainly depends critically on the scat- 
ter of the power-law slope. Our model constrains the 
Compton-thick AGN through mid-IR data which de- 
pends on the SMBH IR/X-ray relationship. The addi- 
tional difference between their work and this study is 
that we included the reflect ed Component for the quasars 
(e.g. iPiconcelli et al.ll2010l ). A summary is given here: 

(1) The unobscured SMBH SED is composed of three 
components, a primary power-law with a cut-off at 
high (> 100 keV) energies, a reflected component and 
iron 6.4 keV emission line. The power-law component 
arises directly from the hot corona above the accretion 
disk. I ts shape is characterized by a photon index of 
r=1.9 (INandrafc Poundlfl99i It" ccvcs fc Turner! [20001: 
IPiconcelli et al.l 120051: iMainieri et al.l 20071) arid a c ut-off 
energy of 300 keV (|Molina et al.ll200alDadindl2008l ). For 
the latter, the uncertainty is still important. This com- 
ponent can be reflected off the far side of the accretion 
disk or dusty torus, hardening the emerging radiation. 
The reflected component depends on the incident spec- 
trum slope and the observed inclination. It is modelled 
for r=1.9 and the relative normalization is fixed to be 1.3 
which is the average for type-1 lines of sight assuming a 
torus half-opening angle of 45°. The iron line is assumed 
to have a Gaussian profile with width of 0.4 keV and 
equivalent width of 280 eV. 

(2) For the obscured Compton-thin AGN, the intrin- 
sic power-law is the same as that for the unobscured 
one but the reflected component normalization is re- 
duced to 0.88 given different viewing angles toward ob- 
scured nuclei. The extinction curve is the photoelec- 
tric absorption at solar abunda nces as characterized by 
iMorrison fc McCammonl (Qj983). In total, we adopts 
log^j=[21.5, 22.5, 23.5] to sample the Compton-thin 
X-ray SEDs. 

(3) For Compton-thick AGN with N m < 10 25 cm" 2 , 
the intrinsic power-law is the same as that for the unob- 
scured one while the normalization of the reflected com- 
ponent is fixed at 0.37. The extinction curve needs to 
consider the non-relativis tic Compton sc attering as mod- 
elled through the code of lYaqoobl (|1997| ). 

(4) For Compton-thick AGN with N m > 10 25 cm" 2 , 
all transmitted photons are down-scattered by Compton 
recoil and only the reflected Component is visible. 

The relative normalization of the X-ray spectrum with 
respect to the IR part is derived through the modi- 
fied version of the tight relationship between the 2-10 
keV and 12.3 urn lu minosities for 22 Seyfert galaxies in 
I Gandhi et al.l (|2009( ). These galaxies are well- resolved 
at 12.3 /xm so that star formation can be removed re- 
liably. Galaxies in that work are a mix of Seyfert 1, 
Seyfert 2 and Compton-thick AGN, while what we need 



here is the relationship only for type-1 AGN, specifically 
PG quasars. At 12.3 /im, the silicate emission feature 
can be important. A s show n in Figure [4j we thus mod- 
ify the iGandhi et a l. (2009) relationship by keeping the 
slope but modifying the normalization to produce the 
median X-ray/12. 3/Ltm ratio of 3 4 PG quasar s whos e X- 
ray and IR data are derived from lCappi et al.l (|2006f ) and 
I Shi et all (|2007l ). respectively. For these PG objects, the 
contribution from star formation to the 12.3 [im emission 
is negligible as evidenced by extreme low aromatic fea- 
ture equivalent width (EWn.3 Mm pAH ;$ 0.05 /im). The 
derived relationship is: 



log 



2-10kcV 

10 43 erg/s 



0.901og 



Ll2.3fj.rn 

10 43 erg/s 



0.45 (15) 



where L2-iokoV an d ^i2.3^m are unobscured 2-10 keV 
and 12.3 /j,m luminosities, respectively. To evaluate the 
effect of this scatter on the model, we explored two addi- 
tional variants of the model with the normalizations 0.2 
dex higher and lower, respectively, where 0.2 dex is the 
3-ct dispersion of the mean ratio of the PG objects. 

3.4.3. HI Column Density Distribution 

The HI column density distribution is key to under- 
standing AGN physics and the cosmic X-ray background. 
Its proper characterization needs unbiased sample se- 
lections and secure measurements of HI column densi- 
ties. In the local universe, optically-selected samples of- 
fer the least bias compared to UV/X-ray (biased against 
heavily-extinguished objects and low-luminosity AGNs) 
and IR (biased ag ainst unobscured o bjects and low- 
luminosity AGNs). iMaiolino fc Riekd (|1995aD have at- 
tempted to compile all Seyfert galaxies at mg < 13.4. 
The Seyfert 2 to Seyfert 1 ratio is found to be as high as 
4:1. This ratio is consistent with that found in the Palo- 
mar nu clear spectrosco pic survey of local galaxies at tub 
< 12.5 (|Ho et al.lll997f) . Both studies dealt carefully with 
host galaxy light dilution and the result should be quite 
reliable. Both studies probed very low-lumin osity AGN 
with b olometric luminosities around 10 9 L©. Hao et al.l 
(2005) instead found 1:1 for Seyfert-2/Seyfert-l ratio 
around Lboi ~ 10 9 L© and decreasing toward higher lu- 
minosities. The discrepancy can be partly due to the 
lu minosity depende nce but also that the SDSS spectra 
in IHao et a l. ( 2005) with large aperture and relatively 
low quality suffer signifi cant host dilution. By bu ilding 
a Seyfert 2 sam p le from | Maiolino fc Riekd (|1995af) and 
IHq et all (|1997D . iRisaliti et alj (|1999f ) measured the HI 
co lumn density distr ibution. Combining this with works 
of ICappi et al.l (|2006f ) for Seyfert 1 galaxies, we charac- 
terized the column density distribution with three bins: 



pNHI _ 



dN 



dlogNm 



0.5/t y p e _i for log; 

= { 0.5(1 -/t ypo -l)(l 

0.5(1 



€ [0,22] 



/ct) for log£%€ [22,2- 
/t vP c-i)/cT for log^e [24,26] 
(16) 



where P Nhi is the probability per logarithm of the HI 
column density, the unobscured AGN fraction /t ypc -i is 
defined as the number of objects with log J^-i < 22 to 
the total, while /ct is defined to be the fraction of CT 
(log^fe > 24) among the obscured AGN (22 < log^fe 
< 26). 
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Beyond the local universe, it is increasingly difficult to 
directly measure the HI column density distribution. Our 
model thus allows both luminosity and redshift evolution 
of /type-i and fcT- 



/typc-l — /o,typc-lX 



1 



xexp( 



-(l0gi?i H R 



(1 + 2)^,type 

9.0)x^ /:typo _ 



(17) 



fcT — /o.CTX- 



1 



(18) 



(1 + Z)^.CT 

xexp(-(logL?f R - 9.0)xA, C t) 

where /o,ty P e-i is fixed to be 0.2 and /o,ct is fixed to 
be 0.5 a t Lttb = 10 9 L^. in order to appro ximate the 
result of IMaiolino fc Riekel (|1995al ). Itfoet all (|1997t ) and 
IRisaliti et al.l ()1999[ ). We also impose a maximum and 
minimum fractions of 95% and 5%, respectively, for each 
column range (type-1, Compton-thin type-2 & Compton- 
thick) . There are four free parameters in total to describe 
the HI column density distribution. 

3.4.4. The UV/optical Part Of Unobscured SED 

The composi t e UV /optical/near-IR SED is taken from 
iRichards et~aH (|2006T ) based on the SDSS type-1 quasars. 
We combined this part with X-ray emission through the 
X-ray/optica l lumi nosity relationship as established by 
ISteffen et al.l (pfWih : 



a ox = -0.1071og( 



L 

t/,2500A 

ergs/s/Hz 



) + 1.739 



(19) 



where a ox =-0.3841og(Z^ ^^/-L^kev) 



3.4.5. The Extinction Curve At UV '/optical /IR Wavelength 

The extinction curve for the UV/optical/IR emission 
is adopted to be that of the Small Magellan Cloud from 
[Pel (fl99l which is a good approx i matio n for reddened 
AGN as discussed in lHopkins et al.l (j2007). To relate the 
amount of extinction in the optical/IR to that in the X- 
ray, we need to assume a gas-to-dust ratio represented by 
Nm/A v in AGN. Studies show that in AGN, the extinc- 
tion in the optical and IR is significantly less than what 
is expected fr om the gas-to- dust r atio in the Milky Way 
or SMC ( e.g. | Risaliti et al.1 [2000t IMaiolino et aTTl200lt 



iShi et all 12006 ). Two possible reasons are that there is 
a large amount of dust-free gas surrounding SMBH or 
that the X-ray and IR emitting sources are spatially off- 
set so that they are not extinguished exactly by the same 
material. We here increased the gas-to-dust ratio by a 
factor of 100 to match the observed average ratio between 
X-ray HI colum n density and silicate absorption feature 
(|Shi et al.H2006h . 

3.5. A Summary Of Model Components 

Our model first assumes total IR LFs that are de- 
scribed by a two power law with four parameters (break 
luminosity and density along with the faint and bright 
end slopes) . The IR LF is assumed to evolve with redshift 
in both break luminosity and density, each of which is a 
polynominal function of redshift with three parameters. 
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Fig. 5. — Examples of SEDs at three given total IR luminosities. 
Note that the model adopts a total of 2340 SEDs at a given total IR 
luminosity and a redshift that covers the range of the SMBH energy 
fraction, the extinction and the scatter in the star-formation SEDs. 
Here it shows SEDs at lower resolution grids over the above three 
parameters for clarity. Different colors are for different extinctions. 

Therefore the total IR LF has a total of 10 parameters. 
The model then introduces the SMBH energy fraction in 
the total IR band to decompose the emission into SMBH 
and star-forming components. At a given IR luminosity, 
the distribution of SMBH energy fractions is assumed to 
be a Gaussian function with both mean and scatter de- 
pending on the luminosity, while the mean is further as- 
sumed to show redshift dependence. There are total five 
free parameters including two for redshift dependence of 
the mean and one for its luminosity dependence, one for 
the scatter at the total IR luminosity of 10 123 L Q and 
one for its luminosity dependence. The mean fraction 
at 10 12 3 Lq is fixed based on local ULIRG objects. For 
each power source (star- forming and SMBH), the model 
then convert the emission in the total IR band to other 
wavelengths using star-forming and AGN SEDs. For the 
former, luminosity dependent empirical SEDs are used. 
For the latter, HI columns are added to extinguish the 
type-1 empirical SEDs. HI columns are divided into 
three ranges representing type-1, Compton-thin type 2 
and Compton thick, which are described by two param- 
eters to give the relative frequencies of occurrence (the 
sum of which needs to be normalized). Each parameter 
has a local value fixed at the SMBH IR luminosity of 10 9 
Lq, and is assumed to follow a simple power law depen- 
dence on each of luminosity and redshift. Therefore the 
HI columns have four parameters in total. The model in 
total has 19 free parameters. 

4. NUMERICAL APPROACH 

Our approach to constraining the model is to fit both 
X-ray and IR survey data, with the basic assumption 
that all the galaxies are experiencing both star formation 
and nuclear accretion activities with relative intensities 
that vary with both redshift and luminosity. The model 
is composed of four components (see §[3] for details): the 
total IR LF, the SMBH energy fraction in the IR, the 
star-formation SED and SMBH SED. In principle, the 
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combination of these four can predict any observed prop- 
erties that are related to galaxy number densities and 
SEDs, such as number counts, redshift distributions, LFs 
etc. 

The model has a total of 19 free parameters. To incor- 
porate the uncertainties on the observational constraints 
that we adopted for the model, we present four vari- 
ants of the model as summarized in Table [T] varying the 
SMBH IR/X-ray ratio and redshift evolution in the star- 
formation SED. The numerical calculation starts with 
the generation of a multi-dimensional array of the flux 
density at a given filter for an object at the given total 
IR luminosity, star-formation color, redshift, SMBH en- 
ergy fraction and HI column density, i.e., /^(filter, Z*tir, 
star-formation-color, z, / BH , Ahi), where the numerical 
grid of each dimension is pre-defined and thus /„ is fixed. 
Each /„ has an associated total probability function (the 
product of four model components): 

P tot = $(L tir ,z)xP bh xP sf " co1oi xP Nhi (20) 

P tot contains all free parameters whose best values will 
be obtained by fitting to the data. The general procedure 
to derive f v and P tot is as follows: 

(a) $(Ltir, z) is a function of Ltir and z. The numer- 
ical grid of Ltir covers the luminosity range from 10 7 Lq 
to 10 14 L & with a resolution of 0.2 dex. The redshift grid 
starts with z=0.001 and increases with Az/z=0.1 until 
z = 1.0 with a total of 63 bins. At z > 1.0, Az is fixed 
to be 0.1 up to z—6. 

(b) The given Ltir is divided into the star-formation 
component and SMBH component 1-xir according 
to the SMBH energy fraction log/ BH . The log/ BH has 
a range from -6.0 to 0.0 with a resolution of 0.2 dex. 
At log/ BH = -6.0, the SMBH contribution at any wave- 
length is negligible. Also, such a low value is needed to 
assure the possibility of negligible SMBH radiation in a 
IR luminous object, i.e., L 2 _iokoV < 10 42 erg/s. 

(c) For the given X^ir, the star-formation SED is as- 
signed according to all possible star-formation color as 

defined by log( {"^as-woiZI) ) • The star-formation color 
varies from Co-3a c to Co+3ct c with a step of 0.5er c where 
Co is the function of I^ir ( see § 13.3. Each star- 
formation color has a probability of p SF - color . With the 
assigned star- formation SED, the (filter) is obtained 
by convolving the SED with the filter curve. Note that 
in the X-ray, the flux is the direct integral of the SED 
between the edges of the band. 

(d) For a given L^f R , the SMBH SED is assigned ac- 
cording to the HI column density. The Am grid is defined 
to be logA ffl /cm" 2 = [0, 21.5, 22.5, 23.5, 24.5, 25.5], each 
of which has a probability of P Nm . The assigned SED is 
then used to derive the L BH (filter). 

(f) The final luminosity in a filter is the sum of 
Lf (filter) and L BH (filter). This is divided by the dis- 
tance modulus to give the final /„ (filter). 

At a given total IR luminosity and redshift, the above 
grids produce 2340 SEDs covering the range of the SMBH 
energy fraction, the extinction and the scatters in the 
star- formation SED, with examples shown in Figure [SJ 
In total, there are 10 million SEDs to cover the IR lumi- 
nosity, redshift, the BH energy fraction, the star forma- 
tion SED and the HI columns. After producing /„ (filter , 



-Ltir, star-formation-color, z, / BH , Ahi), the free pa- 
rameters in P tot are derived by fitting to the observa- 
tions. The fitting scheme is the Levenberg-Marquardt 
least-square techn ique with the IDL program MPFIT.pro 
(|MarkwardUl2009]) . 

5. DATA FOR THE FIT 

Our model fits a total of 617 data points from the lo- 
cal LFs plus two basic types of data from deep surveys, 
namely number counts and redshift distributions. The 
basic information of these data-sets is listed in Table [5J 
The main advantages of our fit strategy compared to 
those fitting directly to LFs at different redshifts include: 
(1) that we isolated the number counts as an independent 
dataset for the fit. These are the most direct result of the 
deep survey observations and contain the least modeling 
uncertainties; (2) that we avoided the uncertainties and 
systematic differences associated with the K-corrections 
for the LFs derived in different studies. 

While a large number of data sets are available in the 
literature, we generally preferred those with good statis- 
tics. We thus included either those based on the com- 
bined fields or individual results of fields with either large 
areas or high sensitivities. We also avoid double fitting 
datasets of the same field but from different authors or 
papers. Overall, the number counts include X-ray data 
at 17-60 kcV, 15-55 keV, 2-10 keV and 0.5-2 keV bands, 
and IR/submm data ranging from 24 (im to 1200 /xm. 
The redshift distributions also cover a range of wave- 
lengths, field areas and limiting flux cuts. For the same 
field, the redshift distributions with different flux cuts 
offer constraints on the model LFs in the different lumi- 
nosity and redshift ranges. The redshift completeness in 
data sets we used is high, > 90% either from spec-z or 
photo-z or both. The error bars include Poisson noise 
and cosmic variances with the latter estimated by two 
point correlation funct ion from direct observatio ns or the 
cosmological model of iTrenti fc Stiavellil (|2008D . Except 
for a few data-sets with high flux cuts, photo-z measure- 
ments always constitute an important fraction especially 
at z > 1. For X-ray sources, photo-z errors can be signif- 
icant due to the faint optical/near- IR counterparts along 
with the power-law featureless shape. We thus added 
additional errors to redshift distributions for the photo-z 
part through the simple bootstrapping method. 

5.1. Number Counts 
5.1.1. X-ray Number Counts 

All s ky surveys in the 17-60 keV by INTE- 
GRAL (iKriyonos et all 12(111 and 15-55 keV by Swift 
(|Aiello et all 120121) probe bright AGN down to 6- 
7xl0 -12 erg/s/cm 2 , whose number counts offer useful 
constraints on the model's behavior in the local uni- 
verse. The Cha ndra deep survey d a ta are used to probe 
high-z universe. [Georgakaki s et al.l ([2008) combined sev- 
eral surveys including CDF-N/S, ECDF-S, EGS, EN1 
and XBOOTES fields, to offer number counts in 0.5-2 
and 2-10 keV with high S/N across a large flux range. 
It reaches 1.0x10" 17 and 1.0x10" 16 erg/s/cm 2 in two 
bands, respectively. Although galactic stars are not re- 
moved in this wor k, their contributions are small (<10%; 
iBauer et al.ll2004D. The additional data in the COSMOS 
field bv lElvis et al.l (|2009) are also included for the fit. 
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Fig. 6. — Fit dataset I - number counts. All symbols are the observations while those in the black are not used for the fitting. The four 
variants of the model are shown with different line styles: solid line for reference variant, dashed line for fast_evol_SED_SF, dot-dashed line 
for low_IR2X_BH and three-dot-dashed line for high_IR2X_BH. The reference model is the one with the minimum \ 2 ■ The fast_evol_SED_SF 
variant assumes strong redshift evolution of the star-forming SED. The low_IR2X_BH and highJR2X_BH variants have X-ray luminosities 
at given IR luminosities for the SMBH radiation 0.2 dex lower and higher, respectively (see Table[T]l. 

5.1.2. Spitzer Number Counts 

We fitted to the Spitzer number counts a t 24, 7 
and 160 ^m measured by iBethermin et al.l (|2010aT) . 
which are generally consistent with previous works 
dPapovich et al.ll2004 iShupe et al.ll2008t iLe Floc'h et al.l 
120091 : iFraver et al.l 120091 ) . The study incorporates data 
from fields including SWIRE, COSMOS and FIDEL. Be- 
sides standard data reduction and photometry measure- 
ments, the incompleteness and Eddington bias are also 
corrected in that work, with further removal of stars. 
The counts reach 35 /xJy, 3.5 mJy and 40 mJy at 24, 
70 and 160 /im, respectively. Their stacking data at two 



The recent 4 Ms CDFS data from lLehmer et al.l (|2012f > 
are further used to probe number counts at the faintest 
end, for which 4-8 keV data along with 0.5-2 keV and 2-10 
keV are all included. An additional lar ge area ASCA sur - 
vey in the 2-10 keV is also included (jUeda et al.l 120051) . 
The small difference (~1%) in the counts between two 
different galactic latitudes of this survey indicates a small 
contribution from Galactic contaminators. The pub- 
lished errors of the number counts in above studies have 
taken into account of Poisson noises and cosmic vari- 
ances. If a fixed Photon index (e.g., 1.4) is used to con- 
vert the count rate to the flux, we added an additional 
uncertainty of 0.1 dex quadratically. 
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Fig. 6. — Continued 

longer wavelengths are not used in the fits but plotted for 
illust ration. At 24 p,m , we ad ded the IRAS 25 ^m data 
from lHacking fc Soiferl (|1991[ ) by assuming a flat SED in 
^/„, which provides constraints in the local universe. 

5.1.3. IRAS Number Counts 

In addition to the IRAS 25 ixm number counts, the 
counts at 60 /im from lLonsdale et al.l (|1990f) are also in- 
cluded to provide constraints at the brightest flux end 
(i.e. local universe). The error is caused by the large 
scale structure of the local universe. We estimated it 
roughly as the difference in the number counts between 
South and North galactic cap. 



Our fit includes the number counts at 250, 350 and 500 
M m of detected sources in the HerMES field as estimated 
by I Oliver et al.l (|2010l ). The depth reaches around 15 
mjy at three wavelengths. Their published uncertainties 
account for Poisson noises, cosmic variances and uncer- 
tainties associated with flux measurements and incom- 
pleteness corrections. Deeper counts obtained through 
P(D) (probability of deflection) analysis of the same field 
(Gl enn eFaIll20Tol) are included only for illu stration. At 
100 a nd 160 ^m, we adopted the result of iBerta et al.l 
(|20T1 from three fields (COSMOS, Lockman Hole and 
GOODS-N). The data reach down to 2.0 mjy at 100 /xm 
and 5 mjy at 160 /im. The estimated errors on the counts 
include Poisson noise and field-to-ficld variations. 



5.1.4. Herschel Number Counts 



5.1.5. SCUBA/ Aztec/Mambo Submm Number Counts 
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Fig. 6. — Continued 

The SCUBA co unts at 850^m ar e colle cted from 
iBorvs et all (|2003T ) and iCoppin et all (|2006l ) with the 
overall flux range from 2.5 to 25 mJy. The Aztec 
nu mber counts at 1 . 1 mm are from two r ecent stud- 
ies, EitiuSeiFin (HoHl) and lScott etaLl (|2012f) . The 
combined counts cover the flux range from about 0.5 to 
10 mJ y. T he Mambo 1.2 mm cou nts are from lGreve et al.l 
(|2004f ) and lLindner et al.l (|201 If ) . 

5.2. Redshift Distributions 

5.2.1. X-ray Redshift Distributions 

For 503 objects in th e Chandra CDF-N 2 Ms catalog 
(| Alexander et al.ll2003D. the redshift measurements are 



com piled from iBarger etaLl (|2005f) . iDonlev et al.l (|2007f ) 
and iTrouille et al.l (|2008T ). In the hard X-ray band, 
we constructed two redshift distributions above limit- 



ing fluxes 3x 10 16 and 3x10 15 erg/s/cm 2 , respectively, 
and in the soft band (0.5-2 keV), above two limiting 
fluxes of 6xl0 -16 and 6xl0~ 17 erg/s/cm 2 . For each 
band, the first flux is near the detection limit and the 
second one is about 10 times brighter. The spectroscopic 
redshift and all redshift (spec-z+photo-z) completeness 
is around 65% and 90%, respectively. 

For 462 source s in the Chandra CDF-S 2Ms catalog 
(ILuo et al.ll2008f), t he red shift measurements are avail- 
able in ILuo et abl (|2010| ). Two redshift distributions 
in the hard-xray are measured above 3.0xl0 -16 and 
3.0xl0 -15 erg/s/cm 2 , respectively. In the soft band, the 
limiting fluxes of 5.0xl0~ 17 and 5.0xl0~ 16 erg/s/cm 2 
are used. The spec-z and all redshift completeness is 
around 55% and 97%, respectively. 

With deep multiple optical/near- IR follow-ups, virtu- 
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ally all 1761 X-ray AGN in the COSMOS field are identi- 
fied (|Civano et alj |2012). The redshift distribution above 
3xl0 -15 erg/s/cm 2 in 2-10 keV and 6xl(r 16 crg/s/cm 2 
in 0.5-2 keV are used for the fit, with spec-z and all red- 
shift completeness reaching 60% and 98%, respectively 

The exten ded CDF-S (ECDF -S) contains a total of 
762 sources (jLehmer et alj|200 5j). The redshift distribu- 
tion in 2-10keV is cons tructed at > 2.0 xlO -15 erg/s/cm 2 
(|Silverman et al.ll2010D . The spec-z and all redshift com- 
pleteness is 64% and 95%, respectively. The XMS offers 
318 X-ray sources at bright fluxes (> x!0 ~ 14 erg/s/cm 2 
) wit h high optical identification (90%) (jBarcons et al.1 
120071) . The redshift distribution is constructed above the 
2-10 keV flux of 5xl0 -14 erg/s/cm 2 with high spectro- 
scopic completeness (87%). 

In the soft band, two additio nal surveys by ROSAT 
are included. iHenrv et al.l (|2006l ) have defined a statisti- 
cal sample above 2xl0~ 14 erg/s/cm 2 from the ROSAT 
north ecliptic pole survey. The spectroscopic redshifts 
are highly complete (>95%, Gioia03). We have measured 
the redshift distribution above 5xl0 -14 erg/s/cm 2 from 
this sample. The second one above 3xl0~ 14 erg/s/cm 2 
is con structed from ROSAT-RIXOS survey (jMason et al.l 
2000) . Note that for these two surveys with high limiting 
fluxes, the galactic stars are identified and removed. 

5.2.2. IR/submm Redshift Distributions 

In the COSMOS field, iLe Floc'h etafl (pool ) have 
measured redshift distributions of 24 /im sources. Here 
we included two measurements with flux limits of 0.08 
and 0.15 mJy, respectively. The one with the highest 
flux limit (0.3 mJy) is not used since the bright sources 
excluded from that study may affect the distribution 
at z < 1. The photo-z completen ess reaches 94%. In 
the GOODS-N. IBarger et all (|2008[ ) carried out a highly 
spec-z complete survey. At /24/im > 0.3 mJy, we mea- 
sured the redshift distribution with spec-z completeness 
of 90%. From the Spitzer 5MUSES data (Helou et al. 
2011), we measured the redshift distribution of 24 fxm 
sources at /24pm > 5 mJy where the spec-z complete- 
ness reaches 98%. At 60/im in the local universe, the 
spec-z distribution of IRAS RBGS sample is included 
dSanders et al.H2003D. 

iChapman et al.1 (|2005l ) carried out a spec-z survey of 
submm sources at 850 /jm and modelled the redshift dis- 
tribution for /sso/im > 5 mJy. We fitted this distribu- 
tion in our study. We further included the redshift dis- 
tribution of s o urces at Aztec 1100 /im as measured by 
ISmolcic et al.l (|2012l ) where the spec-z and redshift com- 
pleteness are roughly 41% and 100%, respectively. 

5.3. Local LFs 

We fitted four lo cal LFs, including l ocal 15-55 keV LF 
(|Aiello et al.ll2012D . l ocal 2-10 keV LF (lUeda et al.ll20H . 
the IRAS 12 ^m LF (|Rush et al.l fl99l. and IRAS 25 /xm 
LF (jShupe et al.lll998h . 

6. DATA FITTING RESULTS 

The fits to the number counts, redshift distribution 
and local LFs are shown in Figure |H Figure [7] and Fig- 
ure [51 respectively, along with the best-fit parameters 
in Table [3] for four variants of the model. The model 
fits 617 data-points in total, with 371, 195 and 51 points 



from number counts, redshift distributions and local LFs, 
respectively. Overall our model reproduces the observa- 
tions reasonably well. The differences between the fits 
and data appear at some data-points but are generally 
smaller than 0.2 dex or within 2-er. Cosmic variances or 
photo-z errors are significant for some of these deviations. 
The reduced x 2 of fo ur variants of the model ranges from 
2.7 to 2.9, with at least 0.75 contri bute d by caveats asso- 
ciated with data themselves (see S 16.41 for details). In the 
following, we go through each version of the model and 
point out the differences between the fit and the data 
and their causes if possible. 

6.1. Fitting Number Counts 

For the reference model (the solid line), the differences 
between the data and the fit include 0.1 dex constant 
over-fits of the 15-55 keV number counts, 0.2-0.3 dex 
under-fits of the 2-10 keV counts at the bright end, the 
0.2 dex under-fits around 0.5-2 keV fluxes from 1.5 to 
4.0xl0 -17 erg/s/cm 2 , 0.1-0.2 dex under-fits around the 
peak 24 /im count in the flux density range of 0.2 to 0.5 
mJy, 0.2 dex under- and over-fits at the 70 /im bright and 
faint end, respectively. Unlike the 15-55 keV counts, the 
17-60 keV ones are well produced by our model so that 
the discrepancy seen in 15-55 keV can be due to subtle 
differences in the X-ray SED around the range of 15 to 
17 keV and 55 to 60 keV. It is yet interestin g to notice 
that, similar to ours, several other models (|Gilli et al.1 
l2007tlTreister et alj|2009t iDraper fc Ballantvnell2010D all 
overproduce the 15-55 keV counts by 0.1-0.2 dex, while 
their adopted SEDs are not exactly the same. The dif- 
ference seen at the bright end of the 2-10 keV counts 
is likely caused by co ntamination of g alaxy clusters in 
the data themselves (|Gilli et al.l 119991 ). In 0.5-2 keV 
counts, the data clearly show a transition flux (3xl0~ 17 
erg/s/cm 2 ) below which normal galaxies other than AGN 
start to dominate, causing the counts to go above the 
simple power law extrapolation from the higher fluxes 
([Bauer et al.ll2004t iLehmer et al.ll2012l ). while the model 
also produces a similar feature but at fluxes ^3 times 
lower. For the 24 /xm peak counts, although the high- 
est observed point i s the result of the SWIR E data near 
the detection limit ( Bethermin et^l2010ah . the differ- 
ences are seen over several adjacent data points. Since 
this peak is caused by aromatic features moving into the 
24 /jm filter around z—1, a possible explanation is that 
z^l star- forming galaxies have slightly increased aro- 
matic feature EWs. Differences seen in 70 /im are also 
noticeable, for which the subtle caveat in star-forming 
SEDs can be the potential cause too. 

For the fast_evol_SED_SF variant (dashed line), ad- 
ditional differences between the fits and data, besides 
those for the reference model, are the significant under- 
estimates of the 70, 100 and 160 /im counts above the flux 
density of the peak count (Figure |6jh)-(j)). The reason 
is that the redshift evolution in the star-formation SED 
results in a colder star-formation color feo^m/ /loo^m 
that decreases the rest-frame 30-100 /im emission. The 
behavior of the other two variants (lowJR2X_BH and 
high_IR2X_BH) in fitting the data are quite similar to 
that of the reference variant. 

6.2. Fitting Redshift Distributions 
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Fig. 7. — Fit dataset II - rcdshift distribution. All symbols are the observations. The four variants of the model are shown with different 
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for high_IR2X_BH (see Table [IJ. 



The redshift distributions are shown in Figure [7] The 
distributions at 2-10 keV are plotted in Figure EJa)-(g). 
The data are from five fields including CDF-N, CDF- 
S, ECDF-S, COSMOS and XMS that complement each 
other by probing from small deep field (0.12 square de- 
gree) to wide shallow area (3 square degree). The lim- 
iting fluxes range from >3xl0 -16 up to > 5xl0~ 14 
erg/s/cm 2 . Compared to the number counts, the obser- 
vational uncertainties here are larger due to the redshift 
measurement errors or cosmic variances. For example, 
for the CDFN and CDFS (Figure (a) and (b), respec- 
tively) with the same field size and limiting X-ray flux, 
we saw obvious differences in the data. In the soft X- 
ray (Figure EJh)-(n)), the data extend from >5xl0~ 17 



to >7xl0 -14 erg/s/cm 2 while the field area ranges from 
0.12 to 80 square degree. Again, cosmic variances and 
photo-z errors play roles in differences among redshift 
distributions with similar flux cuts, emphasizing the lim- 
itations of current available data. The redshift distribu- 
tions in the IR bands are shown in Figure [Tfo)-(u). Al- 
though the redshift measurements in general have high 
completeness, the incompleteness is always redshift de- 
pendent, missing more objects at higher z. This seems to 
be consistent with that our model produces more objects 
than the observations at high z ends in several panels 
(e.g., FigureHm), (n) & (r)). 

All variants of the model perform in similar ways. The 
fit to the 2-10 keV distributions are generally acceptable 



16 



Shi et al. 



2.0 
1.5 
1.0 
0.5 
0.0 
-0.5 
-1.0 



0.2 


1 2 


3 


: (g)2- 


-lOkeV (XMS) 




( 

Li— 1" 


> ^.00E-14erg/s/cm 2 




















spec- 


z: 87% 




L oll-z: 


87% 

■ i 





0.2 



1 



4 5 



N 



cn 
O 



0.2 



0.2 
1 



0.4 
2 3 



0.6 



4 5 



(i)0.5-2keV (CDFN) 

[-1 7erg/s/cm 2 




spec-z: 62% 
oll-z: 86% 



0.2 



0.2 



0.4 
1 2 



0.6 
3 



0.8 
5 



(k)0.5-2keV (CDFN) 

, >^.QCjE- 1 6erq/s/cm 



spec-z: 64% 
oll-z: 87% 

0.2 




4.0 
3.5 
3.0 
2.5 
2.0 
1.5 
1.0 

3.5 
3.0 
2.5 
2.0 
1.5 
1.0 
0.5 

3.5 
3.0 
2.5 
2.0 
1.5 
1.0 
0.5 



(h)0.5-2keV (CDFS) 

f„ > 5.00E-17erg/s/cm 2 




(j)0.5-2keV (CDFS) 

f„ > 5.00E-16erg/s/cm 2 




(l)0.5-2keV (COSMOS) 

f„ > 6.00E-16erg/s/cm 2 




0.4 



0.6 



0.8 



spec-z: 58% 
oll-z: 98% 

0.2 



0.4 



0.6 



log(1.0 + z) 



Fig. 7. — Continued 

where the distribution with the brightest cut as shown in 
Figure [TJg) shows the largest difference between the fits 
and data. In the 0.5-2 keV, two distributions with the 
brightest flux cut (Figure [TJm) and (n)) have worst fits. 
Models over-produce z > 2 sources, potentially caused 
by the fact that sources without redshift measurements 
lie at high z. At 24 /xm, for two lowest flux cuts (Fig- 
ure HJo) & (p)), all variants of the model re-produce well 
two peaks around z=l and z=2 when different aromatic 
features are redshifted into the 24 /im filter. The only dif- 
ference is seen at z < 0.5 where the COSMOS field does 
not probe a significant cosmic volume. For the GOODS- 
N (Figure IHq)), all the fits are systematically lower at z 
< 0.7. Here the cosmic variance or photo-z error should 
also be important; the GOODS-N distribution with a 
higher flux cut but shows higher densities at z < 0.5 



compared to two COSMOS distributions. The match is 
excellent to the highest flux cut at z < 2 (5 mJy in 24 
/xm, Figure Wlj)) where the data quality is high (effec- 
tively 10 squar e degree field siz e and 98% spectroscopic 
completeness, iWu et all I2010D . For this distribution, 
sources with no redshift measurements are not included 
but likely lie at high z according to their photo-z, which 
at least partly causes the difference between the model 
and data at z>2. All the fits to redshifts of local IRAS- 
60/im sources brighter than 5.24 Jy are excellent. All 
the fits to the SCUBA-850 /im and Aztec-1100 fim are 
acceptable but produce higher mean redshifts. 

6.3. Fitting Local LFs 

For the 15-55 keV LF (Figure [8] (a)), all variants of 
the model produce 0.3-0.5 dex higher spatial densities 
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Fig. 7. — Continued 

than observations below the break luminosity. The fits 
to the 2-10 keV LF (Figure (b)) offer acceptable fits 
given the error bars. For the IRAS 12 fxm as shown in 
Figure [8] (c), all variants of the model produce slightly 
flatter bright-end slopes. At the faint end, all variants 
under-estimate the observed LF. The local supercluster 
is known to boost the spatial density of faint sources, 
an effort which is yet not corrected in the 12 am LF 
(|Rush et al.lU99l iFang et al.l[l998t IShupe et al.lll998h . 
For the IRAS 25 /im LFs as shown in Figure HJd) , the 
fits are excellent. 

6.4. Systematics and Model Variants 

All variants of the model, including the reference, 
fast_evol_SED.SF, lowJR2X.BH and high_IR2X.BH as 
listed in TableQ] offer similar performance in fitting num- 



ber counts, redshift distributions and local LFs, with 
the fast_evol_SED_SF variant showing additional signifi- 
cant under-estimates of the 70, 100 and 160 /im number 
counts above the flux density of the peak count. The 
reduced \ 2 ranges from 2.7 to 2.9. The contribution to 
the reduced \ 2 horn systematics associated with data 
is estimated to be at least 0.75 as shown in the follow- 
ing. For the number counts, we assumed negligible data 
systematic errors so that the difference between model 
predictions and data are due to the limitation of the 
model. For the local 12 /j,m LF (Figure HJa)), the ex- 
cess faint sources due to the local super cluster are not 
corrected, causing higher densities than our model's pre- 
dictions. We quantitatively estimate its contribution to 
\ 2 by summing the square of the difference between ob- 
servations and predictions divided by the observed errors 
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Fig. 7. — Continued 

at the 12 /mi luminosity below 10 9 L©. Cosmic variances 
or photo-z errors are evident for redshift distributions 
between different fields with similar flux cuts including 
at 2-10 keV, the CDFN in Figure (a) vs. CDFS in 
Figure [3(b), the ECDFS in Figure [3(c) vs. CDFN in 
Figured] (d), the COSMOS in Figure [7J (e) vs. CDFS 
in Figure [7] (f); at 0.5-2 keV, the CDFS in Figure [7] (h) 
vs. CDFN in Figure (i) , the CDFN in Figure [7j(k) vs. 
COSMOS in Figure [7J (1) . In addition, in 24 fim bands 
at z<0.3, the COSMOS above 0.08 mJy in Figure [7] (o) 
vs. the GOODS-N above 0.3 mJy in Figure (q). Al- 
though this pair has quite different flux cuts, the cosmic 
variance or photo-z errors must play dominant roles be- 
cause the GOODS-N with much higher flux cuts shows 
significantly higher source densities at z < 0.3 compared 
to the COSMOS field. Since the best fit to two Y values 



at fixed X is the mean of two Y values, independent of 
any model assumption, we estimated the minimum con- 
tribution to x 2 for each pair as listed above by assuming 
the best-fit to be the mean of two distributions in the 
pair. The sum of all above contributions normalized by 
the degree of freedom gives 0.75, indicating that the true 
reduced \ 2 that our four variants of the model would 
achieve if data systematics was removed can be better 
than 2.0. 

Comparison of the model variants and their perfor- 
mance in fitting the data does not favor strong star- 
formation SED evolution; on the other hand, direct com- 
parisons with the observed high-z SED color do not pro- 
vide definitive evidence about SED evolution. Figure 
shows such comparisons where we sampled the model 
output to mimic the observational selection, which is cru- 
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Fig. 8. — Fit dataset III — local IRAS and X-ray luminosity functions. All symbols are the observations. The four variants of the model 
are shown with different line styles: solid line for reference variant, dashed line for fast_evol_SED_SF, dot-dashed line for lowJR2X_BH and 
three-dot-dashed line for high_IR2X_BH (see Table [TJ. 

for follow-up observations and thus subject to possible 
contamination from SMBH radio emission. In Fig. [9jc) , 
the tight empirical relationship between i bserved-24^m 
and iobserved-70/im of z=1.3-1.8 star-forming galaxies is 
defined by stacked data and it favors strong evolution in 
the star-formation SED. For the stacked color of z =1.8- 
2.3 galaxies, the variants of no SED evolution performs 
better at the low luminosity end while the variant with 
strong SED evolution does better at the high luminosity 
end. The well-defined relationship between rest-frame 8 
/im and total IR luminosities of the Herschel GOODS 
sample in Fig. Ele) is based on individual detections but 
from different redshifts. Only at the high luminosity 
end, corresponding to z^2, the variant with strong SED 
evolution produce results closer to the observations. In 
Fig- Mi ). the £ bs-24 Aim vs. £ o bs-250 Al m relationship is 
based on 24 /im sources (>0.1 mJy) with Spitzer mid-IR 
spectroscopic measurements at the redshift range of 0.3 
to 2.5. We selected objects with inferred AGN fractions 



rial to eliminate s election biases ([Magnelli et al.1 120091 : 
IChapin et al.ll20lTI ). However, it is impossible to ascer- 
tain that all selection biases have been properly repre- 
sented. Fig. G3a) verifies that the model re-produces the 
local relationship (Equation ITT]) . The difference in the 
median color at a given IR luminosity is around 0.02, 
which is reasonable given the limited resolution (0.037) 
of the numerical calculation. Beyond z=0, we tested 
the evolution of the SED sha pe by comparing to fiy e 
samples, a SMG sample fro m iChapman et al.l (120051) . 
two star-forming samples from Magnelli et al. (2011), the 
Herschel GOODS sample fromlElbaz et al.1 (l201llT and a 
star-forming sample from iKirkpatrick et al.1 (J2012'). For 
the SMG sample, the radio emission in our model is 
derived through the radio-FIR relation (for details, see 
iDale fc Heloul l2002t) . As shown in Fig. [^b), although 
the scenario of no evolution in the star-formation SED 
offers a slightly better match to the 850/Ltm/1.4GHz ra- 
tio of the SMG sample, these SMGs are radio selected 
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(o) IRAS Sample (f MMm >1 .2Jy) ^ (b) SMG Sample 
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Fig. 9. — The comparison in the SED between the model and observations. (Equation [12} • For every observed dataset, we apply similar 

selection functions to our model to define the comparison sample: (a) The local rest-frame IR color (log( jf 6 ^*"" )) vs. IR luminosity of 

IRAS galaxies with feo^m > 1-2 Jy (Chapman ct al. 2003). (b) The observed color between 850 /an and radio 1.4 GHz for the SMG sample 
with /850/jm > 3 mjy, /i.4GHz > 30 /ijy and z > 1 (Chapman ct al. 2005|). (c) The observed 70 fim and 24 fim lumino sity relation of 
star-f orming galaxies at 1.3 < z < 1.8 where open squares are stacking results and filled squares are individual detections l IMaenelli et al.l 
1201111 . Star-forming galaxies in the model are defined by having the fractional SMBH luminosity in the total IR band smaller than 20%. 
(d) The observed 70 /an and 24 /an luminosity relati on of star-for ming galaxies at 1.8 < z < 2.3 where open squares are stacking results 
and filled squares are individual detections (Magnclli ct al. 20 111). (e) The rest-frame IRAC 8 fim and total IR luminosity relation for 
Herschel GOODS sample as defined by fioo um > 1.1 mjy (Elbaz ct al. 2011). (f) The observed 250 fim and 24 fim luminosity relation of 
star-forming galaxies at 0.3 < z < 2.5 (Kirkpatrick ct al. 2012). The four variants of the model are shown with different line styles: solid 
line for reference variant, dashed line for fast_evol_SED_SF, dot-dashed line for low_IR2X_BH and three-dot-dashed line for high_IR2X_BH 
(see Table [TJ. 



smaller than 20% to define the above relationship. The 
250 /an detection rate is quite high (95%) so that the 
relationship should not bias against low 250/im/24/im 
sources. The model variants with no SED evolution 
match better the median of the distribution, compared 
to the variant with SED evolution. At the bright 24 fj,m 
luminosity end that corresponds to 2=2, modelled galax- 
ies with SED evolution are too cold to produce enough 
observed frame 250 /im (i.e., rest- frame 80 /im) radiation. 
We therefore conclude that although the high-z ULIRGs 
are colder than local counterparts as argued in the liter- 
ature, their IR SEDs do not represent exactly those of 



local nor mal star-forming galaxie s across the whole IR 
band (e.g. iKirkpatrick et al.ll20T2h . 

7. THE BASIC OUTPUTS OF THE MODEL 

We discuss here the model's first results in the form of 
its three basic outputs, i.e., the SMBH energy fraction, 
the HI column density and the total IR LFs whose best- 
fit parameters are listed in Table [3] The first of these is 
the unique output compared to previous models for X-ray 
or IR data alone. As a final check of our overall outputs, 
Figure [10] and Figure [TT] show the comparisons between 
data and model's predictions for CIRB and CXB spectra, 
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Fig. 10. — The cosmic IR background light. The red lines 
are model's predictions on the total radiation while black lines 
are just for the emission from SMBH dusty tori: solid line for 
reference variant, dashed line for fast_evol_SED_SF, dot-dashed 
line for lowJR2X_BH and three-dot-dashed line for high_IR2X_BH 
(see Table [TJ. All others are the observational constraints. 
References: Bethe rminlO - IBethermin et al.l d2010al l: Elba z02 - 
lElbaz et all 11200211: Laga c heOO - ILagache fc Pugeti I EOOOh: La- 
gache0 3 — [Lagachc ct al. (2003); Matsuurall - Matsuur a et al 
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Fig. 12. — The median SMBH energy fraction in the total IR 
band as a function of the IR luminosity at four redshifts.The four 
variants of the model are shown with different line styles: solid line 
for reference variant, dashed line for fast_evol_SED_SF, dot-dashed 
line for lowJR2X_BH and three-dot-dashed line for high_IR2X_BH 
(see Table [TJ. 
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Fig. 11. — The cosmic X-ray background light. The red lines are 
model's predictions: solid line for reference variant, dashed line for 
fast_evol_SED_SF, dot-dashed line for lowJR2X_BH and three-dot- 
dashed line for high_IR2X_BH (s ee Table [TJ. All oth ers are the ob- 
servational constraints . jello08 - Aicllo ct al. (2008); Churazov07 - 
IChurazov et al.l (1200711 ; Gend reau95 - IGendreau et alj||1995l); Gru - 
ber99 - IGruber et all JT999I1: M o retti0 9 - IMoretti et al.l 120091); 
Revni vtsev03 -[Rcvnivtscv ct al. (2003|); TurlerlO - ITiirler et al.l 
i201Ch . 

respectively. All variants of the model re-produce the ob- 
served CIRB data within error bars. It is shown that the 
emission from the SMBH dusty torus is around 15% at 
5-10 p and decreases to ward longer wavele ngths, lower 
than Treister et al.l <|2006T ) but comparable to Silv a et all 



(|20(M h All variants of the model reproduce the CXB 
spectrum below the peak energy but are about 10-20% 
lower than observations at higher energies. 

7.1. The SMBH Energy Fraction 

As listed in Table [3] and shown in Figure and Fig- 
ure [T21 all variants of the model indicate that the median 
SMBH energy fraction in the total IR band increases 
with the IR luminosity, which is consistent with vari- 
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Fig. 13. — The redshift evolution of the BH energy fraction 
at Ltir[Lq]=12.3. Different line styles are for different vari- 
ants of the model: solid line for reference variant, dashed line for 
fast_evol_SED_SF, dot-dashed line for lowJR2X_BH and three-dot- 
dashed line for high_IR2X_BH (see Table [TJ. 



ous s t udies in the local universe (e.g. iSanders fc Mirabel 



1996; Gcn zel et al.lll99a lArmus et al. 120071: iDesai et al 
2007tlVeilleux et al.|l2009tlPetric et alBoilD and at high- 
z ( Kartaltepe et al1l2010f l . while the Gaussian width of 
the distribution increases with the decreasing luminos- 
ity. The overall trend of the median SMBH fraction with 
redshift shows a minimum around z~l-2. To test the ro- 
bustness of this redshift trend, we first ran a model with- 
out redshift evolution (kf^ =0, kf^ =0) while all others 
are the same as those in the reference variant. The re- 
sulting best- fit reduced x 2 increased by as much as 1.3, 
which offers strong evidence for the requirement of red- 
shift evolution in the SMBH energy fraction. We then 
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ran a variant of the model that only allows monotonic 
evolution (i.e., /cf^=free, fc™ =0). The increase in the 

best-fit reduced \ 2 is smaller but still significant by 0.4. 

In Figure 1141 we compared the 6.2 /an EWs in our 
model to several observed distributions at different red- 
shifts, which offers a stringent check for our SMBH en- 
ergy fraction measurements. Here the modelled samples 
are defined in similar ways to those observed. All EWs 
are measured through the spline interpolations, with the 
pure star-formation median value around 0.6-0.7 /xm. 
Note that our SF templates do not have scatters in the 
EW, and thus the modelled EWs do not have values in 
bins higher than EW=0.6-0.7 /im and any value below 
EW=0.6 /im is caused by SMBH contributions in our 
model. At logLTiR=12.3 and z=0, the SMBH energy 
fraction in the model is fixed to that of 1-Jy ULIRGs, 
and is allowed to be luminosity and redshift dependent. 
As shown in Figure ITlT a) for ULIRGs, the model pro- 
duces an almost flat trend, r oughly consistent wit h the 
result of GOALS's ULIRGs (IStierwalt et al.ll20ll) that 
additionally show excess objects in the the lowest EW 
bin above a flat trend. For local LIRGs+ULIRGs, the 
distribution has a peak at the SF EW with a tail toward 
the low EW en d, which is roughly consistent with the 
GOAL'S result ()Stierwalt et alJl2012f >. Beyond the local 
universe, we compared th e model's results to those of 
5MUSES (IWu et al.ll2llToh . At 0.05<z<0.3, the higher 
star-formation peak in the model compared to the obser- 
vation is due to zero scatter in our SEE) templates. Over- 
all, the star-formation peak plus a low EW tail as pro- 
duced by the model is consistent with the observations. 
At 0.3<z<0.5, the model reproduces a bimodal EW dis- 
tribution as observed, a SF peak plus a peak at the lowest 
bin with the latter higher than the former. For two higher 
redshift bins (0.5<z<1.0 and 1.0<z<2.0), the model re- 
produces a peak at the lowest bin with a tail toward the 
star-formation end as observed. In Figure H4l we also list 
the fraction of objects with EW < 0.4 for each observed 
or modelled distribution. This quantitative comparison 
further supports the consistency between model's predic- 
tions and observations, with the former is only slightly 
lower than the latter. 
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Fig. 14. — The 6.2 /an EW distributions. The cyan and grey filled 
histograms represent the observed detections and upperlimits, re- 
specti vely. The local LIRGs and ULIRGs are from Sticrwalt et 
(2012), while those at higher redshifts are from 5MUSES (Wu et 
12013) . The open histograms are the model's predictions. Differ- 
ent line styles are for different variants of the model: solid line for 
the reference variant, dashed line for the fast_evol_SED_SF variant, 
dot-dashed line for the low_IR2X_BH variant and three-dot-dashed 
line for high_IR2X_BH (see Table [j). The listed percenta ges are 
the fractions of objects with EW < 0.4 for each distribution ("obs" 
stands for the observed distribution while "mod" is for the mod- 
elled distributions). 




7.2. HI Column Density Distribution And 
Compton-thick AGN Fraction 

The model characterizes the HI column density dis- 
tribution through three bins, the type-1 (logTVui/cm 2 < 
22), Compton-thin type-2 (22 < logA ffl /cm 2 < 24) and 
Compton thick AGN (24 < logiVHi/cm 2 < 26). Figure 
[TBI shows the evolution of the type-1 and Compton-thick 
AGN fractions in all AGN as a function of the SMBH 
luminosity and redshift. All four variants of the model 
require that the type-1 fraction increases with the BH lu- 
minosity but decreases with redshift, and the CT fraction 
has the opposite trend. 

The trend of the type-1 fraction with the SMBH lumi- 
nosity is in general consistent with previous studies (e.g . 
Ueda et al.l 120031 : lHasingerl 120041 . lLa Franca et al.l 120051: 



Barger et alJl2005b iTreister et alJl2006t iBallantyne et al.l 
20061: iGilli et al.ll2007t lHasingerl 120081) . Although previ- 



ous CXB models can explain the X-ray survey data ei- 
ther with the redshift dependence of the type- 1 fraction 
(|La Franca et alj 120051 IBallantyne et al.l I2006T ) or with- 



FlG. 15. — The redshift evolution of typc-1 AGN fraction 
(left panel) and Compton-thick AGN fraction in all AGN pop- 
ulations. Different line styles are for different variants of the 
model: solid line for the reference variant, dashed line for the 
fast_evol_SED_SF variant, dot-dashed line for the low_IR2X_BH 
variant and three-dot-dashed line for high_IR2X_BH (see Table [TJ. 
Different colors are for different X-ray luminosities: blue for 
log(L 2 _iokev[erg/s])=42.5, cyan for log(L 2 _i 0feeV [erg/s])=43.5, 
green for log(L 2 _iokev[ er g/ s D = 44.5, orange for 

log(L 2 _iokev[erg/s])=45.5, red for log(L 2 _i k c v [erg/s])=46.5. 

out it (|Ueda et al.H2003HTreister fc UrrvlllX)! IGilli et al.l 
120071 ). the requirement of such a redshift evolution is 
highly preferred by our model that is constructed to ex- 
plain both the X-ray and IR data, consis tent with recen t 
studies of a large X-ray AGN sample (|Hasingerl 12008) . 
Without redshift evolution (/3 z ,t yP c-i =0 & Pz,gt =0), 
the best-fit reduced \ 2 rises by 0.5, which is essentially 
caused by the increased brightness of high-z AGN at both 
X-ray and 24 /zm wavelengths. In the pure CXB model 
that only deals with the X-ray data, the increase in the 
apparent brightness of individual objects can be comple- 
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mented by the decrease in the spatial number density 
in order to keep the number of objects above a given 
flux density unchanged. This is why the CXB model is 
not sensitive to redshift evolution of the type-1 fraction. 
However, the change of the AGN number density in our 
model is additionally constrained by the IR data. The 
redshift evolution of the break luminosity and density of 
the total IR LF is largely fixed through IR survey data. 
Although the change in the SMBH energy fraction can 
vary the AGN number densities, it also modifies the num- 
ber density of the star-formation-dominated sources es- 
pecially the far-IR/submm sources that lie at the bright 
end of the IR LF and thus are sensitive to the SMBH 
energy fraction. Without redshift evolution of the HI 
column density, the best-fit result show that the redshift 
distributions at z > 1 of 24 /im sources brighter 5 mJy 
and 0.5-2 keV sources brighter than 10 -14 erg/s/cm 2 are 
over-produced by a factor of 10, and the redshift distribu- 
tion of submm sources systematically offsets by a factor 
of 3 along the Y-axis. 

Our model requires a large number of the Compton- 
thick AGN. At < 10 keV, Compton-thick AGN are ex- 
tremely faint. The traditional CXB model first con- 
strains the Compton-thin AGN through X-ray data at < 
10 keV and then add the Compton-thick AGN to match 
the CXB spectrum above 10 keV especially t he 20-30 keV 
peak fe.g. lUeda et~aLl[2003t [Gilli et al.ll2007h . In this pa- 
per, we take advantage of the capability of the mid-IR 
data in probing these deeply buried AGN thanks to sig- 
nificantly lower apparent gas-to-dust ratio. A test run 
of the model with zero Compton-thick AGN abundance 
shows that the reduced x 2 would worsen by as much as 
1.0. The main reason is that the relative fraction be- 
tween Compton-thick and Compton-thin type-2 AGN as 
a function of the flux density is different at 24 /im from 
that at 0.5-2 and 2-10 keV, so that it is able to break 
the degeneracy between Compton-thick and Compton- 
thin type-2 AGN seen in pure CXB models. As shown in 
Figure [TBI the Compton-thick AGN are only important 
over the faintest 1-2 orders of magnitude flux ranges in 
X-ray bands, while Compton-thin type-2 AGN are im- 
portant over a much larger flux ranges. The resulting 
Compton-thick to Compton-thin type-2 AGN fraction in- 
creases with the decreasing flux density at both bands. 
Such similar behaviors at both soft and hard bands im- 
ply that if the model is forced to replace Compton-thick 
objects with Compton-thin type-2 AGN in one band, the 
result in the other band would not significantly violate 
the observation. However, at 24 /im, the Compton-thick 
contribution is noticeable over a flux range as large as 
four orders of magnitude, which is quite similar to that 
of the Compton-thin type-2 AGN. The relative ratio be- 
tween the two populations is not a strong function of 
the flux density, with a value roughly around 1:3. This 
means, if the model is forced to replace Compton-thick 
AGN with Compton-thin AGN at 24 /im, the result in 
X-ray bands would significantly over-produce objects at 
the high flux end but undcr-producc at the faint flux end. 
This argument was verified by running a test model with- 
out a Compton-thick AGN population and the expected 
behavior was observed in the model output. 

7.3. The Total IR LF 
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Fig. 16. — The percentage energy contributions to the number 
counts at 2-10 keV, 0.5-2 keV and 24 /im, from different types of 
radiations in the reference model. Cyan for the pure star-formation 
emission, yellow for the pure SMBH emission of type-1, orange for 
the pure SMBH emission of Compton-thin type-2 and red for the 
pure SMBH emission of Compton-thick AGN. 

As shown in Figure [T71 the predicted total IR LFs are 
consistent with the observations over redshift and lumi- 
nosity ranges where the data are available. For compar- 
ison, we converted the observed LFs at given filters to 
the total IR LFs using the reference variant. We first 
modelled the number density ^modci (filter, z, vL v ) at 
the same filter, redshift bin (z±Az) and luminosity bin 
(vL v ±AvL v ) as the observed ones (<E> bs (filter, z, vL u )). 
The vL v is then corrected to the total IR luminosity 
with luminosity corrections based on the probability- 
weighted average of all model SEDs within z±Az and 
vL v ±AvL v . This < I ) modci/ < I ) obs is taken as the difference 
of the observed and modelled LFs in the total IR band. 
For the local LFs, the IRAS data offer the largest dy- 
namical range and are thus the best dataset to test the 
validity of the model's LF. The 12 and 25 /im LFs are 
used to fit the model and the result is discussed in § 16.31 
Here we further compared our model to t he total IR LF 
based on IRAS 60 /im selected sources ([Sanders et alj 
2003), and found a good match between the two (Fig- 
ure [T7]). Consistency is excellent with other works of sig- 
nificantly smaller fields that cover the whole IR/submm 
wavelength range. Above z > 0.1 up to z ~ 3, cur- 
rent deep surveys are able to probe objects up to one 
order of magnitude below the break lum inosity as shown 
in Figure [171 At short wavelengths rtWu et all 120111: 
Rodiehiero et al.ll201Ct iRuiopakarn et al.ll2010l) T no sig- 
nificant difference is found between the model and data 
up to z^2.5. Some deviations at the faint end are most 
likely caused by incompleteness corrections particularl y 
associated with 1/V max method ()Page fe; Carreral feOOO). 
A t rest-frame 250 /im , the model is co nsistent wi t h that 
of lEales et all (flOlOh but the result of iDve et all pOlOl ) 
is systematically higher at 0.1< z <0. The total IR LF of 
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Fig. 17. — The luminosity function in the total IR (8-1000/^m) at different redshifts. Different line styles are for different variants of 
the model: solid line for the reference variant, dashed line for the fast_evol_SED_SF variant, dot-dashed line for the l ow_IR2X_BH varian t 
and thr ee-dot-dashed lin e for high.I R2X.F3H (see Tab le [J. References : Caputi07 -ICaputi et al.l (120071) : Dunn eOO -IDunne et al.l 112000) : 
DyelO - IDve et all lt2010l1 : E aleslO - lEales et all II2010T): Rodigh ierolO - IRodighiero et al.l (120101 ) : Ruiopakarn lO -rRuiopakarn et all fcoid) : 
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IWu et al.l72011l) . 

and then gradually reaches the peak at z=4, while the 
break density shows a gradual increase until z=0.4 and 
then monotonically decreases with redshift. To test the 
redshift out to which luminosity or density evolution is 
required, we have carried out several test runs of the 
model by modifying the evolution formula to force no 
evolution above redshifts of z—0.5, 1.0, 1.5, 2.0, 2.5 & 
3.0. If adopting Ax 2 /d.o.f=l as an indicator of signifi- 
cant worse fit, we found that the no-evolution scenario 
can be ruled out below z ~ 2.5. Above z~3, the fit result 
could not differentiate evolution vs. no evolution, which 
is not surprising since the current constraint on the z > 
3 infrared and X-ray populations is still of low statistical 
significance as shown in Figure [7J To evaluate the de- 
generacy between the luminosity and density evolution, 
we carried out test runs which halted density or lumi- 
nosity evolution at selected redshifts. The result shows 



ICaputi et al.l ((2007) is systematically above our modelled 
LF; their work may over-estim ate the total IR lumi nosity 
due to K-corrections (also see IMaenelli et al1[20rl . 

Our model hypothesizes no evolution in either the faint 
or bright end slopes of the total IR LF. As shown in Fig- 
ure [171 the fixed bright end slope matches the available 
data up to z~3. For the faint end slope, we did not see 
deviations from the data up to z^l above that the avail- 
able data are quite limited to justify. With fixed faint 
and bright end slopes, the total IR LF as predicted by 
our model requires evolution in both luminosity and den- 
sity. Although the detail relies on the evolution formula 
adopted, we found in general that the break luminos- 
ity (logL*) increases with redshift but the break den- 
sity (log^*) decreases with redshift. As shown in Fig- 
ure [TSl for the evolution formula given in Equations [2] 
and [3l the break luminosity increases rapidly until z—2 
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Fig. 18. — Redshift evolution of four parameters describing the total IR LF including the break luminosity, break density, faint end 
and high end slopes. Different line styles are for different variants of the model: solid line for the reference variant, dashed line for the 
fast_evol_SED_SF variant, dot-dashed line for the lowJR2X_BH variant and three-dot-dashed line for high_IR2X_BH (see Table [TJ. 



that the data are not able to break the degeneracy at z 
> 1.5, which is consistent with Figure [171 showing that 
above z ~ 1.5 the observed LF does not extend signifi- 
cantly below the break luminosity However, no matter 
which evolution formula is adopted, the robust result is 
that the break luminosity always increases with redshift, 
while the break density overall decreases with redshift. 
Our result of the positive luminosity evolution is consis- 
tent with different literature studies while the negative 
trend in the density is not always favored. We note our 
result of negative density e volut io n is inconsistent with 
studie s in iLe Floc'h et al.l (12005D. [Perez-Gonzalez et al.l 
(2005) and iRodighiero et al.1 (120101) . but consistent wi th 
ones in lCaputi et al.l (|2007j ) and iMagnelli et al.1 (|2011[ ). 

8. CONCLUSIONS 

In this paper, we have constructed a joint model of 
the X-ray and IR extragalactic background to account 
simultaneously for both X-ray and IR survey data. The 



main conclusions are the following: 

(1) The model successfully fits the deep survey data 
with the scenario that individual galaxies are experienc- 
ing both star formation and SMBH accretion activities. 
The model with 19 free parameters generates number 
counts, redshift distributions and local luminosity func- 
tions consistent with 617 data-points over six decades 
of electromagnetic frequencies. The best-reduced x 2 
reaches 2.7-2.9, with at least 0.75 contributed by cos- 
mic variance, photo-z errors or limitations of the data 
sets. 

(2) The unique output of this model, compared to pre- 
vious ones for either pure X-ray or IR data, is to con- 
strain the SMBH energy fraction in the total IR band. 
The best-fit requires the SMBH energy fraction to in- 
crease with the IR luminosity but decrease with the red- 
shift back to z^l.5. The derived trend is statistically 
significant as shown by the significantly increased best- 
fit x 2 /d.o.f., i.e., Ax 2 /d.o.f > 1.0, for best fits where the 
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SMBH fraction has either no luminosity or no redshift de- 
pendence. An independent test comparing the model's 
predictions to the observed aromatic feature equivalent 
width distributions agrees with this result. 

(3) The second important output of the model is the 
distribution of the HI column density obscuring the 
AGN, especially the Compton-thick fraction. A run of 
the model with no Compton-thick AGN shows an in- 
crease in the x 2 /d.o.f. by ~ 1.0. The best-fit parame- 
ters indicate that the Compton-thick AGN fraction de- 
creases with the SMBH luminosity but increases with the 
redshift, while the type-1 AGN fraction has the inverse 
trend. 

(4) The third output is the total IR luminosity func- 
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TABLE 1 

A Summary of Different Variants of The Model 



Models 
[1] 


description 
[2] 


Pc 

[3] 


Lo § L f H 10kcV VS - 

[4] 






reference (solid line) 


reference-model 


o.o : 


| n „- L 2-10kcV _ 


0.901og 




0.45 


IO{ = 10 4a crg/s 


' 10 43 crg/s 


fast_evol_SED_SF (dashed line) 


strong evoluton of SF SED 


2.0 


l„_- L 2-10koV _ 


0.901og 


, L 12 .3u, m 


0.45 


l °S 1043 erg /s 


■ 10 43 erg/s 


lowJR2X_BH (dot-dashed line) 


lower BH IR/X-ray ratio 


o.o : 


1 -E-2-lOkoV _ 


0.901og 


, L 12 .3nm 


0.25 


l °S 10 43 erg/s 


■ 10 43 eig/s 


high_IR2X_BH (three-dot-dashed line) 


higher BH IR/X-ray ratio 


o.o : 


| n „- L 2-10koV _ 


0.901og 


il2.3 M m 


0.65 


lo g iQ4a crg / s 


' 10 43 erg/s 



Note. — Col.(l): The names of different variants of the model; Col. (2): the description to the specific variant of the 
model; Col. (3): the p c value f Equation 1 12|| used to quantify the star-forming SED redshift evolution. Col. (4): the SMBH 
X-ray vs. IR luminosity relationship used to connect the AGN X-ray SED part with the IR part. 

Compared to the reference model, the fast_cvol_SED_SF variant assumes strong redshift evolution of the star-forming SED. 
The low_IR2X_BH and high_IR2X_BH variants have X-ray luminosities at given IR luminosities for the SMBH radiation 0.2 
dex lower and higher, respectively. Sec § |3.3.2| § 13.4,21 and § 16.41 The line style for each variant of the model is used for all 
figures in this paper. 
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TABLE 2 

Data Sets To Derive The Best-Fit Parameters 









differential number counts (371 data 


points) 


band 


ref 


Fields 

All ~1„. 


Area^] Depth 



15-55 keV 
4-8 keV 
2-10 keV 



Aicl lo et al. (2012) 
Lehmer et al. (2012) 



0.5-2 keV 



Gcorgakakis et al. (2008 



Ueda et al. (2005) 
Elvis et al. (2009) 
Lehmer et al. (2012) 
Georgakakis et al. (2008) 



Spitzer-24/^m 



Elvis et al. (2009) 
Lehmer et al. (2012) 
Bcthcrmin et al. (2010a) 



IRAS-25/im 
Spitzer-70/im 



Hacking fc S oifcr (1991) 
Bcthcrmin ct al. (2010a) 



Spitzer-160^tm Betherm in_et_aL Q2010a) 



Herschel-lOO/im Berta et al. (2011 



Herschel-160/xm F3erta_et i ^L i (2011^ 



Herschel-250Atm Oliver et al. (2010) 



Herschel-350Atm Oliver et al. (2010 



Herschel-500Atm Oliver et al. (2010 



SCUBA-850/im 
Aztec-1100/jm 




All-sky All-sky 

CDFS 0.13 

CDF-N 0.11 

CDF-S 0.06 

ECDF-S 0.25 

EGS 0.63 

EN1 1.47 

XBOOTES 9.24 

AMSS 278 

COSMOS 0.5 

CDFS 0.13 

CDF-N 0.11 

CDF-S 0.06 

ECDF-S 0.25 

EGS 0.63 

EN1 1.47 

XBOOTES 9.24 

COSMOS 0.5 

CDFS 0.13 

FIDEL ECDF-S 0.23 

FIDEL EGS 0.41 

COSMOS 2.73 

SWIRE LH 10.04 

SWIRE EN1 9.98 

SWIRE EN2 5.36 

SWIRE ESI 7.45 

SWIRE CDFS 8.42 

SWIRE XMM 8.93 

All-sky All-sky 

FIDEL ECDF-S 0.19 

COSMOS 2.41 

SWIRE LH 11.88 

SWIRE EN1 9.98 

SWIRE EN2 5.34 

SWIRE ESI 7.43 

SWIRE CDFS 8.28 

FIDEL EGS 0.38 

COSMOS 2.58 

SWIRE LH 11.10 

SWIRE EN1 9.30 

SWIRE EN2 4.98 

SWIRE ESI 6.71 

SWIRE CDFS 7.87 

GOODS-N 0.083 

GOODS-S 0.083 

LH 0.18 

COSMOS 2.04 

GOODS-N 0.083 

GOODS-S 0.083 

LH 0.18 

COSMOS 2.04 

A2218 0.0225 

FLS 5.81 

Lockman-North 0.34 

Lockman-SWIRE 13.20 

GOODS-N 0.25 

A2218 0.0225 

FLS 5.81 

Lockman-North 0.34 

Lockman-SWIRE 13.20 

GOODS-N 0.25 

A2218 0.0225 

FLS 5.81 

Lockman-North 0.34 

Lockman-SWIRE 13.20 

GOODS-N 0.25 

HDF-N 0.046 

SHADES 0.20 
GOODS-N,LH,GOODS-S,ADF-S,SXDF 1.6 



6xl0~ 12 erg/s/cm 2 [5<r] 
4.6x10" 17 erg/s/cm [P dct >0.004] 
10" 16 erg/s/cm 2 [P<4xl0" 6 ] 
2xl0 -16 erg/s/cm 2 [P<4xl0" 6 ] 
3xl0 -16 erg/s/cm 2 [P<4xl0- 6 ] 
3xl0" 16 erg/s/cm 2 [P<4xl0- 6 ] 
lO" 15 erg/s/cm 2 [P<4xl0" 6 ] 
6xl0" 16 erg/s/cm 2 [P<4xl0~ 6 ] 
10~ 13 erg/s/cm 2 [5a] 
7.3X10" 16 erg/s/cm 2 [P<2xl0~ 5 ] 
5.1X10- 18 erg/s/cm 2 [P dct >0.004] 
10~ 17 erg/s/cm 2 [P<4x 10" " 



2x10" 
3x10" 



erg 
erg 



;/s/cm 2 [P<4xl0" 
;/s/cm 2 [P<4xl0" 



3xl0" 17 erg/s/cm 2 [P<4xl0" 6 ] 
10" 16 erg/s/cm 2 [P<4xl0" 6 ] 
6xl0 -17 erg/s/cm 2 [P<4xl0" 6 ] 
1.9X10" 16 erg/s/cm 2 [P<2xl0 -5 ] 
5.1 xlO- 18 erg/s/cm 2 [P det >0.004] 
60 /ijy [80% complete 
76 /ijy [80% complete 
96 /ijy [80% complete 
282 ^jjy [80% complete 
261 AiJy [80% complete 
267 /tJy [80% complete 
411 /tjy [80% complete 
281 /iJy [80% complete 
351 fijy [80% complete 
300 mjy [5<r] 
4.6 mjy [80% complete] 
7.9 mjy [80% complete] 
25.4 mjy [80% complete 
14.7 mjy [80% complete 
26.0 mjy [80% complete 
36.4 mjy [80% complete 
24.7 mjy [80% complete 

45 mjy [80% complete 

46 mjy [80% complete 
92 mjy [80% complete 
94 mjy [80% complete 
90 mjy [80% complete^ 
130 mjy [80% complete] 
88 mjy [80% complete] 



3.0 mjy 
1.2 mjy 

3.6 mjy 
9.0 mjy 

5.7 mjy 

2.4 mjy 

7.5 mjy 
10.2 mjy 

13.8 mjy 
17.5 mjy 
13.7 mjy 
25.7 mjy 
12.0 mjy 

16.0 mjy 

18.9 mjy 
16.5 mjy 
27.5 mjy 

13.7 mjy 

15.1 mjy 
21.4 mjy 
16.0 mjy 
33.4 mjy 

12.8 mjy 
~ 3mjy 



2mJy [la] 
0.4-1.7 mjy [la] 



3a 
3a 
'3a 
3a 
'3a 
3a 
'3a 
[3 a] 

50% complete 
50% complete 
50% complete 
50% complete 
50% complete 
50% complete 
50% complete 
50% complete 
50% complete 
50% complete 
50% complete 
50% complete 
50% complete 
50% complete 
50% complete 
ia] 
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TABLE 2 — Continued 



Mambo-1200/jm 



Hatsukadc ct al 
Greve et al. (2004 
Lindner et aL72011) 



ADF-S 
EN2, LH 
LH-N 



0.25 
0.10 
0.16 



~ 0.5mJy [la] 
0.6 mjy [la] 
0.75 mjy [la] 



Bands 
2-10 keV 



0.5-2 keV 



Spitzer-24/^m 



IRAS-60Mm 
SCUBA-850Mm 
Aztec- 1100/jm 



Fields 

CDF-N 

COSMOS 

CDF-S 

ECDF-S 

XMS 

CDF-N 

COSMOS 

CDF-S 

ROSAT-NEP 

ROSAT-RIXOS 

COSMOS 

GOODS-S 

5MUSES 

All Sky 

7 fields 

COSMOS 



rcdshift distributions (195 data points) 



all-z 
"87%- 

98% 
97% 
95% 
87% 
86% 

98% 

97% 

97% 

93% 

100% 

90% 

98% 

100% 

75% 

100% 



FT 



spec-z 
61% 

60% 
51% 
64% 
87% 
63% 

59% 

50% 

97% 

93% 

0% 

90% 

98% 

100% 

75% 

41% 



Bargcr ct al. (2005), Donley ct al. 
Trouille et al. (2008) 
Civano et al. ('2012) 
Luo et al. TSbloT " 
Silverman et al (2010 ) 
Barcons et al. (2007) 
Bareer et aT"72005) . Donley et al. (2007) 
Trouille et al. (2008) 
Civano et al. (20121 
Luo et al. (ToioT " 
Henry et al. (2006) 
Mason et al. (2Wb) 
Le Floc'h et al. (2009) 
Barger et al. (2008") 
Wu et al. {ToibT ~ 
Sanders et a!"T2003) 
Chapman et al. (2005) 
Smolcic et al. (20TTT 



local luminosity functions (51 data points) 



Bands 
15-55 keV 
2-10 keV 
IRAS-12/im 
IRAS-25Mm 



ref 

Aicllo ct al. (2012) 
Ueda et al. (20TTT 
Rush et al. (1993) 
Shupe et al. (1998) 
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TABLE 3 
The Best-Fit Parameters 



Parameter 


reference 


fast_evol_SED_SF 


low.IR2X.BH 


high_IR2X_BH 


[1] 




[2] 


[3] 


[4] 


[5] 


The total IR LF 




1 


10.61±0.04 


10.67±0.04 




1_U . U .\J*± 






0.35±0.40 


1.65±0.37 


1.65±0.36 


0.81±0.39 


k 2 ',i 




15.33±1.63 


6.97±1.50 


9.70±1.46 


13.00±1.57 


k 3,l 




-15.19±1.76 


-7.64±1.63 


-9.32±1.60 


-12.95±1.70 






-2.67±0.06 


-2.71±0.06 


_o fi4-i-n nfi 

. VJrt_l_ U.UU 


_o fifi+n 06 


K l,d 


3.68±0.57 


2.96±0.57 


9 ^7-1-0 
z.o i nzv.Qo 


9 sn-i-n 


K 2,d 




-16.25±2.24 


-10.90±2.25 


1 1 70-1-9 
-11. / UltZ.Uo 


19 1 9-1-9 91 
-1Z.1ZIZ.Z1 






9.39±2.36 


3.61±2.42 


4.77±2.12 


5 33+2 34 


71 




0.39±0.03 


0.40±0.03 


U. OOIU.UO 


u . ouiu . uo 


72 




2.49±0.04 


2.65±0.05 


2.52±0.05 


2.52±0.04 


BH Energy Fraction In The Total IR Band 










0.64±0.03 


0.75±0.05 


0.55±0.02 


0.63±0.02 


"l,d 




-3.33±0.27 


-3.97±0.32 


-3.29±0.22 


-3.15±0.08 


1.BH 

K 2,d 

P/BH 




3.05±0.44 


6.02±0.52 


2.89±0.40 


2.87±0.11 




0.74±0.08 


0.78±0.08 


0.67±0.07 


1.14±0.01 


Po- 




-0.19±0.05 


-0.13±0.05 


-0.09±0.05 


-0.46±0.01 




HI Column Density 






/^2,typc— 1 




1.23±0.18 


0.79±0.17 


0.91±0.17 


0.68±0.13 


Pi, type— 1 




-0.66±0.06 


-0.48±0.06 


-0.48±0.07 


-0.42±0.05 


/3 Z ,CT 




-0.79±0.09 


-1.00±0.09 


-1.36±0.13 


-0.81±0.08 


A.CT 




0.44±0.08 


0.41±0.07 


0.73±0.12 


0.21±0.05 


d.o.f. 




598 


598 


598 


598 


x 2 




1675.7 


1752.8 


1702.6 


1608.1 



Note. - — 



The best-fit parameters are listed for four variants of the mode 1 (see TabkflJ. 



